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ABSTRACT 


\ 


An efficient and user-oriented method has been construct J for calculating 
flow in and about complex inlet configurations. Efficiency is attained by: the 
use of a panel method, a technique of superposition for obtaining solutions at 
any inlet operating condition, and employment of an advanced matrix- iteration 
technique for solving large full systems of equations, including the nonlinear 
equations for the Kutta condition. User concerns are addressed by the provi- 
sion of several novel graphical output options that, taken together, y»“ld a 
more complete conprehension of the flowfield than had been possible previously. 
Examples of these features are presented for some complicated configurations. 
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PRINCIPAL NOTATION 




A Matrix of induced normal velocities at the control points. Also 

flat panel area. 

b Subscripts 32, 41. Intercept on n axis of panel side. 

3 Dipole derivative along an N-line. Subscripts F and S denote values 

on first and second N-lines of a panel, respectively. Subscript K 
denotes value associated with K-th lifting strip (Fig. 1). 

c Geometric constant for a panel designed to minimize dipole d-> on- 

tinuity along a lifting strip. 

( X ) ( y ) 

C K Geometric constants expressing source derivative on a panel in 

terms of values of source density on surrounding panels. 

d Length of a side of a panel; used with subscripts to denote a par- 

ticular side. 

h Geometric constant denoting arc length along an N-line from trailing 

edge to n-axis of a panel. Subscripts F and S denote first and 
second N-lines, respectively. 

i.j Integer subscripts denoting panel ramfeer. 

Unit vectors along axis of Cartesian coordinate system. Subscript e 
denotes panel coordinate system. 

L Number of lifting strips. 

L(total) Total arc length of an N-l ine from trailing edge to leading edge and 
back again. 

m Subscripts 32, 41. Slope of a panel side. 

n L» n n» n c Components of n in panel coordinates. 

Unit normal vector to a panel. 

N Number of panels. 

N-line Curve along which input points are distributed. On lifting portions 
it defines a wing section (Fig. 1). 

P,Q,R Derivatives of panel shape at origin of panel coordinates 

r Magnitude of r 

r Vector between the two points in space. 

S Surface area of curved panel. Used with subscripts 32, 41 denotes 

cosine of slope angle of a panel side. 


v 


t Maximum dimension of a panel. 

U,V sent at ion 0n *^ C00rdinates used in P^araetric cubic surface repre- 

V Velocity. Used with various subscripts and superscripts. 

VlJ panel lty ** *’ th contro1 P oint due t0 uni 't source density on j-th 

Vi Final combined velocity at i-th control point. 

1 Tiftij5 y st^ip' th COntr<>1 P ° int due t0 un1t dipcle derivative on K-th 

w Width of a panel or lifting strip. 

Cartesian coordinates. 

Coordinates of a point of a panel in its own coordinate system. 

M ?** P** 1 that defines the panel 

ZrtJX S pl2T?S^^f* er1,rtI » «* y *n<rte fcr1»rt1«s at the 

0 th^ori^in^f^SILi paneT - t Subscripts x and y denote derivatives at 

orlg" , " t ** r 5ubscr, ' ,ts «•"•** «'* 

“ Vector vortic.ty strength on a panel. 


vi 




1.0 INTRODUCTION 


The present work consists of the construction of a computer program for analyz- 
ing flow in and about very complicated three-dimensional inlet configurations. 
The basic calculational technique is a panel method, whose choice is dictated 
by considerations of numerical efficiency and geometric generality. Such a 
method can calculate flow about virtual!* any passive configuration in a rou- 
tine fashion. An inlet, however, is active in the sense that it contains 
complicated internal machinery that controls the amount of fluid that enters. 
For the present purpose, it is not necessary or even desirable to analyze this 
machinery in detail. Instead, its effect is lumped into a single parameter, 
namely the mass flow through the inlet. This situation is seen most clearly 
in the static case where the inlet is at rest In an otherwise undisturbed 
fluid. Here the only fluid motion Is that due to ingestion of fluid by the 
inlet. The difference between a method applicable to inlets and a method for 
passive bodies is that the former must contain a calculational device for con- 
trolling mass flow through the inlet. As will be seen, this requi re m en t is 
equivalent to a calculational device for generating a static solution. 

Reference 1 describes a method for calculating flow about sinple three- 
dimensional inlets by means of a first-order panel method. Reference 2 pre- 
sents a procedure applicable to inlets having auxiliary inlets which uses a 
higher-order panel method. The present program generalizes this last to the 
case of inlets having leading-edge slats. Thus this method must account for 
lifting effects, while that of Ref. 2 did not. This represents a considerable 
difference. Not only must the code contain formulas giving the effects of 
bound and trailing vorticity, but the iterative matrix solution must be altered 
radically to include the nonlinear Kutta condition. Moreover, tne user must 
now bear the responsibility of specifying the location of the trailing vortex 
wake. As far as the code is concerned, the designation "leading-edge slat" 
may be interpreted very generally as a lifting portion of the configuration. 
Thus the program can consider an inlet mounted on a lifting wing-pylon 
configuration. 

It should be pointed out that the rather elaborate scheme outlined here is made 
necessary by the requirement to simultaneously analyze flow outside and inside 
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the inlet. If only the exterior flow were of interest, the inlet entrance 
could be represented by panels on which an inflow velocity distribution ,s 

specified. 

This inlet program is thus the latest in a series of panel-method programs, 
each of which depends on the previous ones. Two previous inlet programs. 

Refs. 1 and 2, have already been mentioned. The present code is based on the 
higher-order lifting panel method of Ref. 3, which in turn is based on the 
higher-order nonlifting panel method of Ref. 4, and both are ex ensio 
first-order lifting panel method of Ref. 5. Recent documents, notably Refs. 2 
and 3, refer liberally to the earlier documents. Thus no one document contains 

.11 the necessary formulas and logic of the ensile of prog™* tha ‘ 

developed up to this time. It was decided to remedy this 

the preset rep*t c-plete. All the p-el-eth* f— 1» 

ures that pertain to the present code are contained herein. The resulting 

reference Manual is soaarftat lengthy, but it Is complete. 
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2.0 THE HIGHER- ORDER PANEL METHOD 


\ 


2.1 General Description 

The panel method here Is a higher-order source method. In a source method the 
unknown value of source density on each panel is adjusted to satisfy the 
normal-velocity boundary condition at the panel control points, while dipole 
and/or vorticity is treated as an auxiliary singularity that is used to produce 
the lifting effects. In contrast to a first-order panel method (Ref. 5 . which 
uses a constant source density on flat panels, the higher-order method e . 
accounts for source-derivative and surface-curvature effects. As exaqiles wi 
show, inclusion of these effects can be quite l^ortant for internal flows. 

The panel Mthod of Ref. 3 is the only current method that accounts for local 
surface curvature despite analysis that Indicates that a method is not truly 
higher order if these effects are ignored (see Section 2.2 below). 


A three-dimensional lifting flow is characterised by the P«senaof a trailing 
vortex sheet that issues fro. the trailing edge of the lifting portion of the 
configuration, e.g. a win, or slat (see Figure 1). If -or. than one lifting 
device is present, such a sheet Issues from every trailing edge. The location 
of the trailing vortex sheet is not known a priori. In the present panel 
method the location is simply Input by the user based on his experience and 
physical intuition. In the case of a wing-fuselage, this issue is unl-portant 
because the wake is relatively far fro. all parts of the wing and is weak near 
the body, in the inlet case, however, the wake of the slat is ingested an 
probably lies close to the inner surface of the Inlet. The wake location c 
be important, and its estimation might be somewhat difficult, particularly at 

angle of attack. 


A key issue in any lifting procedure is the method of applying the Kutta con- 
dition. Oddly most panel-method publications virtually ignore this question, 
presumably because the authors consider it unimportant. In fact, the Kutta 
condition is of first importance, because it determines the circulation dis- 
tribution that drives the whole lifting flow. Classically, the Kutta conditw 
is stated as the avoidance of an infinite velocity at the trailing edge. 

Clearly such a condition cannot be enforced nimerically, and some other 
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criterion, which is a consequence of the Kutta condition, must be Invoked. 
Physically, it .s absolutely necessary that upper- and lower-surface pressures 
on the wing must approach a ccimnon value at the trailing edge. This is the 
form of the Kutta condition employed in the present panel method. Oddly no 
other panel method uses the equal -pres sure Kutta condition. Instead, other 
derivative conditions are used. The reason for this is probably that the 
equal-pressure condition is nonlinear, while the alternate conditions are 
linear, which simplifies the numerical procedure. The price of linearity is 
high, however, because calculated results can be seriously in error. Figure 2 
shows calculated results from Ref. 6. It can be seen that only the present 
panel method gives equal upper- and lower-surface trai ling-edge pressures for 
this case. The other methods give a pressure mismatch of up to half of free- 
stream dynamic pressure. 

toother i^wrtant feature of the p re s en t panel method Is its use of an Itera- 
tive matrix solution, to inlet having a centerbody, an amllfary Inlet, and 
one or more slats Is a very complicated configuration, which if mounted on a 
wing-pylon becomes much more complicated. Thus the cases to which the present 
method will be applied tend to have large panel nutoers. For such cases an 
iterative solution of the linear equations that express the boundary condition 
is an order of magnitude faster computationally than a direr*, elimination 
solution. As is well known, the computational effort for an Iterative solution 
is proportional to the square of the number of equations, i.e. to the panel 
nuafcer, while that of the direct solution is proportional to the cube. Thus 
the advantage of the former becomes greater as the panel nimfcer increases and 
as the computation tine becomes a more important factor. For the iterative 
solution to realize this advantage, it must converge reliably In a relatively 
small number of iterations, say 10-20. As reported In Ref. 2, an accelerated 
block Gauss-Sledel achieved this efficiency in the nonlifting panel method. 

The addition of lift produces a major complication in that the values of bound 
vorticity must be Included among the unknowns and the set of equations must 
contain those expressing the Kutta conditions at various span locations along 
the trailing edges. As mentioned above, these last are nonlinear. The stand- 
ard "Newton-type" method for solving sets of nonlinear equations consists of 
successive linearizations followed by Iteration. That Is, In each Iteration a 
set of linear equations is solved, but the coefficient matrix changes with 
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each iteration. This alternative typo of iteration was incorporated into the 
basic block Gauss-Siedel scheme, and after a restructuring of the acceleration 
procedure, a reliable method was obtained. Hie number of iterations required 
for convergence is increased 30% compared to the nonlifting case. 

2.2 Consistency Analysis 

One of the distinguishing features of panel methods is that the velocity due 
to a panel singularity is computed analytically, as contrasted with other 
so-called "boundary-element" methods where numerical quadratures are employed. 
This feature is rendered necessary by the requirement faced by a practical 
panel method that panel dimensions are often larger than characteristic physi- 
cal dimensions of the boundary. For example, in the midchord mid-semi-span 
region of a wing the spanwise dimension of panels is normally an order of mag- 
nitude larger than the local wing thickness (Fig. 1). The integrals over a 
panel giving the potentials and/or velocities due to various singularity dis- 
tributions can be integrated analytically only if the panel Is flat. Other 
panel methods have assumed flat panels from the beginning and ther hyp othesize d 
singularity distributions that are two-variable polynomials of various degrees, 
usually ranging from constant to quadratic. Investigators who use the poly- 
nomials of higher degree tend to label their methods "hi*er-order" and thus 
to imply that the increase of accuracy with panel nurter of such a method is 
more rapid than that of a "lower order" method. Such as assertion has been 
proven false by direct comparisons of calculated results (Refs. 6 and 7). The 
reason is singly that successive refinement of an integrand (the singularity 
distribution) without simultaneous refinement of the integration region (the 
panel geometry) cannot lead to improved results, because the factors neglected 
are more important than the additional factors included. 

An alternative approach is to expand the effect of a general panel about its 
tangent panel. All relevant quantities can be expanded in Taylor's series 
about the tangency point and the Integral expressed as a sequence of terms each 
of "higher order" in panel dimension than the preceding terms. Since all 
integrals are over the flat tangent panel, they all can be evaluated analyt- 
ically. But the expanded terms contain derivatives of both the singularity 
strength and the body shape. The analysis is done in detail in Appendix A for 
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the case of source singularity, it is shown there that a term of a given order 
contains derivatives of the surface shape that are one degree higher than the 
highest singularity derivative. Thus consistent combinations are: flat-panel/ 
constant source, paraboloidal panel/linear source, cubic panel /quadratic 
source, etc. The same is true for vorticity, which reflects the fact that 
vorticity effects can be expressed in terms of source effects (see Appendix H). 
Dipole effects have a more complicated expansion due to the fact that a dipole 
distribution on a panel is equivalent to a vorticity distribution equal to its 
gradient (and thus a polynomial one degree lower) plus a concentrated vortex 
filament around the edges of the panel. The present approach uses panel vor- 
ticity and adds the appropriate edge vortices. 

2.3 Development of the Panels from Input Points 

As in all panel Methods, the body surface and m i c e are input to the 
by specifying the coordinates of a number of points on the surface. These are 
associated in groups of four to form the quadrilateral surface panels (Fig. 1). 
This may be done in a variety of ways. In the present program the end result 
is a trapezoidal tangent panel (Fig. 3) and various geo m e t r ic quantities (about 
60) associated with it. This is a much smaller number of geometric quantities 
than many other panel methods require. The order of the input points is along 
certain curves called N-lines (Fig. 1). on lifting portions of the body the 
first and last points on the N-line are at the trailing edge. The set of 
panels formed from the points lying on two consecutive N-lines on a lifting 
portion is denoted a lifting strip of panels. 


The initial step in generating the panel consists of using a "canned" routine 
for fitting surfaces by parametric bicubic splines. This is applied to each 
panel individually to generate the panel coordinate system, the coordinates of 
the four comers of the panel in this system, and the three second derivatives 
of the surface at the origin of panel coordinates. The panel coordinate origin 
is the point of tangency to the surface of the tangent panel and is also the 
control point where the normal-velocity boundary condition is applied. The 
procedure for doing this is described in Appendix B. 
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As shown In Fig, 3, the corner point coordinates are £ K » n K » * * 2 » 3 » 4 

The lengths of the parallel sides are: 

<J 12 > dp * d 43 ’ 1 S s5 r t 4 U- 

The width of the panel is 


The slanting sides are straight lines with equations of the form 

e; = mn + b 


(2.3.3) 


where 


h ’h 

*32 * w * 

«l-«4 

*41 = w 

Kf\Z ~ ^3 


b 32 " w * 

d 41 w 

The maximum diagonal of the panel is 



[ / U 2 - S 4 ) Z + (n 2 - n 4 ) 2 


t = max 




(2.3.4) 


(2.3.5) 


Further define 


S 


32 




S 41 



(2.3.6) 


The lengths of the slanting sides are 




(2.3.7) 


Also needed for lifting panels are the total arc lengths along the H-lines from 
th e trailing edge up to the n-axls of the panel In question. These are 
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hp = [ dp - h s = Z d s - r >4 


(2.3.8) 


\ 


where the sums are over the previous panels of the lifting strip. 

Finally the normalized moments of the area of the panel are required. The 
method of calculating these is in Appendix C. 

2.4 Velocities Induced by the Source Distribution on the Panels 

Formulas for the velocity induced by an individual panel at a point in space 
are obtained from the expressions developed in Appendix A by differentiating 
and then performing the indicated integrations over the flat projected panel. 
Different procedures are called for depending on the distance of the point in 
question from the panel. For nearby points, the expressions of Appendix A are 
integrated exactly. This procedure is rather lengthy and details are omitted. 
Only the final formulas for this “near field" are presented (Appendix D) and 
these are the key to the present panel method. It is assumed that the point 
in question has been transformed into the panel coordinate system and afl near 
field formulas are given in terms of this coordinate system. For points 
further from the panel, the integrals of Appendix A are evaluated by a classic 
multipole expansion. The orders of the expansions are selected to be at least 
as high as the terms in question. This is a relatively sinple procedure ana- 
lytically, and the resulting "intermediate field" formulas require much less 
computing time than the near-field formulas (Appendix E). This computation 
also is carried out in panel coordinates. The velocities calculated by the 
near-field and intermediate-field formulas must be transformed into the refer- 
ence coordinate system. For points even further away, a "far-field" approxi- 
mation is used (Appendix F). This is obtained simply by retaining only the 
first terms in the multipole expansions. However, the far-field formulas have 
been put in vector form, and thus they can be evaluated directly In the 
so-called reference coordinate system in which the body is input. This elim- 
inates the need for transformations and further reduces computing time. Some 
of the quantities in the near-field formulas lose numerical significance for 
certain ranges of values of the parameters, e.g. control point near the exten- 
sion of a side or effect of a very long thin panel on adjacent control points. 
Most of these problems are due to the short word length used by IBM computers 
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and do not arise for CDC equipment. A variety of special formulas based on 
power series expansions have been developed for use in the troublesome situa- 

tions. These are collected in Appendix G. 

The source potential is given in Eqs. (A.31) through (A.35) of Appendix A The 

velocity induced by the panel is obtained by taking the negative gradient 

obtain 


, . t<0>o * !PD ,P1 ♦ 2Qt (Q) ♦ rt ,R1 - » (1x, -’x ♦ ?(1y V 


(2.4.1) 


where each individual velocity is the negative gradient of the corresponding 
potential . 

2.5 The Source Derivative Ter; Assenfcly of the Matrix of Influence, 
Coefficients 


2.5.1 The Numerical Differential Procedure. Geometric Constants 

AS stated in Section 2.4, the basic source velocity formula (2.4.1) contains 
coefficients o_ and o , which are the derivatives of the source 
with respect to the £»V. coordinate directions. It is not intend hat 
these be additional unknowns. Instead, they are expressed in terns of the 
unknown values of source density at the control points of the surrounding 
panels Thus, ultimately the values of source density at the control points 
of the panels are the only unknowns. The source-derivative procedure is 
slightly different for the first and last panels of a strip and for pane s of 
the first and last strips. However, the modifications are quite straicfit- 
forward In the initial discussion it is assumed that the panel on which 
source derivatives are being evaluated (the panel in question) has adjacent 
panels on all four sides as shown in Fig. 4. For the purposes of the present 
discussion only, the control points of the adjacent panels are "inhered 1 = 0, 

1 2 3, 4 as shown in Fig. 4, where 0 denotes the element in question. Dies 

control points are transformed into the coordinate system of the panel in 
question. Let the % and n coordinates of these control points be L 0K , ig*. 

K - 0 1 4 and the values of source density at the control points be o K . 

K = o’, l’...,4 (evidently L 00 - igg ' 0) • The differe " tiat1 °" pr ° CeSS < ” < '"' eSS ' S 
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the source derivatives on the panel in question in terms of the unknown values 
of source density on the adjacent panels in the form 


V p{*>„ 

°x = Jo * K 


M 


= l c 

y K=0 


(y) c 

K °K 


(2.5.1) 


where H represents the number of adjacent panels. It is 4 for interior panels, 
3 for panels on the edge of a section, and 2 for panels in a corner of a 

section. 


The essentials of the numerical process are that it calculates one-dimensional 
derivatives in the u and v directions of the parametric-cubic coordinate system 
of Appendix A and then expresses the derivatives with respect to the panel 
coordinates in terms of these. 

For each panel, calculate the geometric quantities 


u = At/a, v = w/a. 

For the purpose of one-dimensional differentiation, define the coordinates 

X 1 = * 7 < 3 0 + 3 1 ^* x 2 = 7 <*0 + 3 2 ^ ^ 

t 3 - \ («o ♦ a 3 >. *4 = "7 (a 0 + a 4 } 

where the subscripts are panel designations of Fig. 4. Then centered 3-point 
differentiation, which is appropriate in the interior of a section, gives 

X 2 x l * *2 X 1 

°x = x 1 (x 2 - °1 " °o ' x 2 (x 2 - x^ °2 

\ c. • 0 • 

t 3 *3 + *4 *4 

°t = t 4 (t 3 r t^T °4 - -Tfc °o ' t 3 (t 3 - t 4 ) 3 


a = / (At) 2 + 


(2.5.2) 


a = 


(dp ♦ d $ ) 
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The first of Eqs. (2.5.4) already gives the first of Eqs. (2.5.1). Thus 


. . x. + x, 

r(x) = _ _1 2 

0 x^x 2 


c (x) = X 2 

1 x 1 1 X 2 * x i ) 


C 


(x) - , 

2 x 2 (x 2 


X 1 



r (x) 

L 3 


r (x) 

l 4 


= 0 


(2.5.5) 


By analogy the second of Eqs. (2.5.4) gives t-derivative coefficients 


„(t) *3 At) _ t 3 + *4 r (t) _ _ t4 

C 4 ' t 4 (tj - t 4 T ’ C 0 * * * 


t 3 t 4 * ~3 yVH T 


By the chain rule 


cj 1 * * C 2 tJ = 0 


0£ = 0 X U + OyV 


Thus 


°y * v (ff t ■ °x u) * I Y " uC K X,)o K 


(2.5.6) 


(2.5.7) 

(2.5.8) 


and finally 

c[ y) = ^ (C,^ - uc£ x) ), K = 0 , 1 , ..., 4 (2.5.9) 

For panels on the edge of sections, the centered 3-point formulas (2.5.4) must 
be replaced by 2-point one-sided formulas in the direction (or directions) 
where a third value does not exist. 


2.5.2 Logic of the Assembly Procedure 

In the first-order method the velocity induced by a panel depends only on the 
source density on that panel and thus the "influence coefficients" for that 
panel are calculated solely from that panel’s geometry. The essentially new 
feature of the source derivative procedure is that the velocity induced by a 
panel depends on the value of source density at the control point of that panel 
and also on the values of source density at the control points of adjacent 
elements. Similarly, the velocities induced by adjacent elements depend on the 
source density on the panel in question. Thus the "influence coefficients" for 
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a panel depend not only on the geo*try of that panel but also on the geometry 
of adjacent panels and the assembly of the influence coefficient matrix is more 

complicated. 

let the panels be numbered consecutively in the order they have been formed. 
Thus, reference is made to the i-th panel and to the j-th panel where both i 
and j range from 1 to N. Another way of stating the essentially new feature 
above is that a distinction must be made between the effect of the j-th E52L. 
and the effect of the j-th value of source den|ity, whereas these two < effects 
are identical in the first-order method, let t*„ be the velocity Indue* at 
the 1-th control point by the j-th panel and be the velocity induced at 
that point by the j-th value of source density. Them n the notation of 
Section 2.4 and the present section, 

. [♦«» * fO’lp a 2l“»Q * * c‘ X ¥ ,X> ♦ ^ ]o„ 

■ ^ 


? [£<*)?««» a 


(2.5.10) 


notice that subscripts i and j are omitted on the right side of Eq. * Z ’*^®* 
for si»11city. In the overall nurtering scheme, o 0 in (2.5.10) is oj o,, 
o, and o 4 have subscripts near j. All velocities in (2.5.10) depend onb- 
on tte geometry of the j-th panel. The curvatures P, Q, R and the coeffici 
ct*) c<» depend on the surrounding panels, but once calculated they can 

be associated with the j-th panel only. 


;onsider now the i-th row of the matrix *4> ich expresses the ®* fe ** s ° f 
the various values of source density at the i-th control point. 

bracketed tern in (2.5.10) is an effect of and is added to the *** 

, r mtr f A mK in the suMMtion of (2*5.10) represent 

location of the row. The four terns in tne swwt.u _ 

effects of other values of o and must be added to other locations. Referring 
to Fig. 4, It can be seen that the panels nurtured 1 and 2 are on the same 
strip as the panel in question and thus represent effects of the preceding and 
succeeding values of o. In particular, value 1 is associated with 03., and 
value 2 with and the relevant terms of Eq. (2.5.10) are added to those oca- 
tions. Panels 3 and 4, however, are on adjacent strips. Suppose there ore 
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„ aneU on „ ch Strip. Then value 3 is associated with 0j . E and valued .1th 

l and the relevant terns of Eq. (2.5.10) are added to these locates 
j+E* 

the row. 

2.6 Vorticity Influence of a Panel 

2.6.1 Panel Vorticity and the Underlying Dipole Distribution 

the effect of a vorticity distribution on a panel cannot be expressed in terms 
of a potential, as mentioned above. The velocity Induced by the panel of Fig. 

8 at a point (x,y,z) is 


u x r 


-d$ 


( 2 . 6 . 1 ) 


( 2 . 6 . 2 ) 


* f J „ 

m s r ~ 

Mhere « is the vector vorticity strength and 

f * (x - E)t + (y - n)j ♦ Cz “ 

The distance r has its usual neanh* «<h1<* *1“ equal to \f\. and the 

integral Is over the true panel. To insure that the 

usual vorticity conservation theorens over the panel, . _ 

express S in terms of equivalent dipole distribution p. As show in Ref. , 

the relation is 

- - . (2.6.3) 

w = -n x grady 

were X is the unit normal vector, .hose comments nj, V are given by 
Eqs (A.8) through (A.10). In Appendix A, T, j, k «ere used as uni ,ec or 
il the axes oTthe panel coordinate system, because no other coordinate 
svTten entered the discussion. Here, to avoid any possible confusion, these 

unit vectors .ill be Witten %, l e , £ e (A**™ 11 * B > t0 SBecify they 
indeed unit vectors of the panel system. To be compatible .1th the source 
£Hty a"d pane, -geometry expansions, u is taken to be a quadratic function 

of panel coordinates E and n in 

(2.6.4) 


JL 2 

y = y 0 + v x S * v/» + v xx r + + Hyy* 
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so that the components of u> vary linearly over the panel. Furthermore, since 
a two-term expansion of the source potential Is all that appears feasible, only 
a two-term expansion of Eq. (2.6.1) is required. Evidently 


& * “x ♦ “"xx 5 ♦ V> 
In * M y + 2 tM xy 5 + ’Vy' 1 * 

are two-terra expansions of the components of grady. 


( 2 . 6 . 5 ) 


Then Eq. (2.6.3) gives 


W 


ln C ' 


(n. 


*L)1 + 
3? ,J e 


(n n \ Jn^e 


( 2 . 6 . 6 ) 


Two- term expansions of the ? e and 3e co^onents of (2.6.6) contain zeroth 
and first-order terms. Since the leading term of the component is linear, 
only that one term is required. The two-term expulsion of Eq. (2.6.6) is 

£ - Im + 2(i*^ + y yy n )J^e * L V 2Cw xx c + v xy n)Jj e + L 'Wx + <2sV*e 

(2.6.7) 


Non using Eq. (2.6.7), a taw-term expansion of Eq. (2.6.1) may be carried out, 
and the resulting velocity put in terms of source influences. This development 
1$ curried out in Appendix H. 


The dipole strength is required to vary linearly over the H- lines bounding the 
panel. In particular 


on n * ru 


( 2 . 6 . 8 ) 


u = Bp(C + h F ) 
y * B $ (E * \) o" n = 

when these are applied to Eq. (2.6.4), it turns out that y must have the form 


y * g Un + hpn - - n^ip ♦ cw(n - njXn - tv| )3®p 

- J Un ♦ h $ n " + cw(n - njMn - n^JBj 


(2.6.9) 
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Notice that y is expressed as the sum of two dipole distributions: one 

multiplying B p that is zero for n = and one multiplying B s that is zero for 
11 =^. Thus, one dipole distribution is associated with each N-line. The 
logic of the calculation keeps these two separate until a later stage of the 
calculation. The constants B p and B $ are essentially bound vorticity strengths 
that are determined by the Kutta condition. It is the distribution multiplying 
each that is important at this stage, so in effect the B's are set equal to 
unity. The constants in the underlying dipole distributions are 


y- derivative 

First N-line 

Second N-line 


n l 

^3 


w 

w 

"y 

hp 

- cti^ + Hj) 

he 

~\T* ctn l + V 

“xx 

0 

0 (2.6.10) 


1 

1 

M xy 

7w 

" « 

Pyy 

c 

-c 


The foregoing are on-body formulas. For wake panels set 


v 


X 



0 


ty = L p (total), h $ = L $ (total) 


( 2 . 6 . 11 ) 


where L (total) is the total arc length of an N-line from trailing edge to 
trailing edge. Equation (2.6.11) reflects the fact that the underlying dipole 
strength is constant along N-lines in the wake. 


All constants in Eq. (2.6.10) are known except c. It is determined to make the 
dipole strength as continuous as possible from one panel to the next along a 
lifting strip. Clearly nothing enforces continuity if there is a physical gap 
between the panels, so c is determined assuming that adjacent panels on a 
lifting strip share a coanon side. This seems the best that can be done. 
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Consider the dipole strength along the "top" side of the panel between the 
points (C 3 , hj) and (^, n 2 ) (Fig. 3). It is obtained by setting K - m 32 n + 
b 32 in Eq. (2.6.9). The result is 

v>( 32 ) = y( linear) + (Bp - B«j){cw^ + wm^} [•£ (^ “ D] (2.6.12) 

In the square bracket s denotes arc length along the side and L the total 
length of the side (L = d 32 in the notation of Section 2.3). The function 
y( linear) is a linear function that varies from the value of y at the point 
(S 3 , nj) to the value of y at the point (s^ t^). On the adjacent element, the 
"bottom" side that lies between the points (S 4 * and (Sj» nj) is the one 
that lies along the side discussed above. The dipole strength along this side 
is 

v(41) « y( linear) ♦ (Bp - BjMcw 2 + wm 41 > [f (£ - 1)1 (2.6.13) 

Ignoring any small gaps b e t wee n elements, the quantities y( linear), s, and L 
are identical in Eqs. (2.6.10) and (2.6.11), as are Bp and B^ The only 
quantities that are different are those in the curly brackets. Here c and w 
correspond to different elements, while the slopes m^ and » 4] correspond 
to different sides of different elements. Thus continuity between panels i 
and i + 1 of a strip is obtained if 


„(i) 


[c( i >w^> + «4j ) ]=w 1+1) [c^' +1 U i+1) 


4 -»♦!)] 

+ *41 J 


(2.6.14) 


where w is panel width (usually the same for all panels of a strip), n 32 is 
the slope of the upper panel edge and m 4 i the slo^e^of the lower panel edge, 


Eq. (2.6.14) is solved for successive values of c 

:< 1 > =0 


beginning with 


(2.6.15) 


and proceeding over all on-body panels of the strip. The choice, Eq. (2.6.15), 
is arbitrary and expresses the fact that Eq. (2.6.14) has a nonunique solution. 
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2.6.2 Edge Vortices 


The fundamental development of Ref. 5 shows that the velocity induced at a 
point in space by a dipole distribution p over a panel is identical to the 
sum of the velocity induced by a vorticity distribution w, as given by Eq. 
(2.6.3), and the velocity due to a concentrated line vortex around the 
perimeter of the panel whose variable strength equals the local value of 
dipole strength. While the velocity due to the dipole distribution is 
inherently a potential flow (zero curl), neither of the other two velocities 
are; only their sum is potential. Using u> in the form of Eq. (2.6.3) satis- 
fies the vorticity conservation theorems over the surface of the panel but not 
at its edges. Thus to the vorticity velocity of Appendix H must be added the 
effects of line vertices on the edges of the panel with strength equal to the 
local value of the underlying dipole distribution. Since an actual body obvi- 
ously does not have line vortices in its surface, in the absence of numerical 
approximation the edge vortices of adjacent panels would cancel exactly. Thus 
it might be hoped that the panel edge vortices could be ignored away from 
physical edges such as wing tips. It turns out that this Is true for spanwlse 
panel edges but not for streamwise panel edges. That is, referring to Fig. 1, 
the edge vortices of adjacent panels on the same lifting strip cancel to a good 
approximation (especially when the continuity algorithm of Section 2.6.1 is 
employed) and thus nay be ignored. However, the edge vortices that lie along 
an N-line in general do not cancel with those of panels of the adjacent lifting 
strip to a degree that justifies their omission. 

There are several ways of accounting for the effect of the edge vortex, all of 
which are theoretically equivalent to some order of accuracy. The approach 
used here is the analogy of that used throughout the higher-order development. 

A vortex lying along the edge of a curved panel is projected into the tangent 
plane. The relevant formulas are developed in Appendices I and J. 

2.6.3 The Trailing Vortex Wake 

The wake is input to the program by specifying points along N-lines just as 
for on-body points. The option exists of making the last panel on each wake 
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strip semi-infinite. In many cases, such as a clean wing, the location of the 
wake has very little effect on the solution. In such cases, the wake may be 
taken as semi-infinite right from the trailing edge, and no wake points need 
be specified. This optional wake may have the direction of either the 
trailing-edge bisector or the x-axis. 

Wake panels have vorticity but no source density. However, because of the way 
in which vorticity effects are calculated in the present program, essentially 
the same induced- velocity formulas must be evaluated as for on-body panels. 

Of course, no boundary conditions are applied on wake panels, and their pres- 
ence does not affect the order of the matrix of the linear equations for the 
source density. 


For finite wake panels, the basic influence formulas are unchanged, but the 
constants defining the underling dipole distribution and the edge vortex form- 
ulas are modified as described In the previous sections. Also modified are the 
values of c that improve dipole continuity between panels of a lifting strip 
(Section 2.6.1). Let superscript (1) denote quantities associated with the 
first on-body element of a lifting strip and superscript u denote quantities 
associated with the last on-body element of the strip. Similarly, the super- 
scripts wl, w2, etc. denote the first wake element, second wake element, etc. 
of the same lifting strip. The important value of c is c 1 , i.e., the one 
for the first wake element. It is confuted from 


c (w1) = 


w( u >[w< u >c (u > + 


Ju)-I _ 
*32 J 
,(wl ) -i2' 


,(l) [w OUl) J I) 


+ m, 


'41 




[w lwu ] 


(2.6.15) 


where the quantities w, m 32 , n 41 have their usual moaning (usually c (1) = 0). 
Values o* c for the remaining wake elements are obtained from 


c (wl) Cw (wl) 3 2 . c (w2) Cw (w2) ] 2 . c (w3) [w (w3) ]2 = ... 


(2.6.16) 


In most cases of interest, the trailing vortex wake extends to infinity. To 
facilitate accounting for this condition, provision has been made for consid- 
ering the last element of the wake to be semi-infinite. A finite element of 
the sort shown in Fig. 3 is formed at the end of the wake, including all the 
geometric quantities of Section 2.3. The induced velocity calculation for this 
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ife 


■ 


i 


a 


element is performed using the origin of coordinates appropriate to the finite 
element, but the formulas used to calculate induced velocities are appropriate 
for the semi-infinite element. Naturally, all points in space are in the "near 
field" with respect to a semi-infinite element, so it is the formulas of 
Appendix D that apply. These formulas are modified by setting 


**32 


= 0 


(2.6.17) 


E 2 ♦ 


^3 00 


This yields inwediately 


(2.6.18) 


otj, P r y r V Y 4 unchanged 

«3 - a 2 = -1, *3 - P 2 * y 3 * y 3 * 0 

The log functions. Eg. (D.3), and their derivatives, Eq. (D.6) are replaced by 


« unchanged, all derivatives unchanged 
1^32) - o, all derivatives equal zero 


(2.6.19) 



-L tl21 ♦ L 1341 

r, - (x - t 4 ) 
* l0 9 r, - lx -qT 

3L (W 

“4 - 1 

a i - ' 

ax 

~ r 4 - (x - E 4 f 

3x r-j - (x - Ej) 

8L (34) 

*4 

6 1 

ay 

' r 4 - (x - E 4 T 

ay r 1 - (x - Ej) 

9 l (34) 

y 4 

a (,2) . T i 

az 

- r 4 - (x - 

az r, - (x - EjJ 


( 2 . 6 . 20 ) 


( 2 . 6 . 21 ) 


The inverse tangent functions, Eqs. (D.4), and their derivatives, Eqs. (D.5), 
are replaced by 
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T< 32 > - tan- 1 ( y ; ' ) 

k = 3 or 2 

(2.6.22) 

t| 1 4 ^ = unchanged, 

3T< 32 > 

^ o 

k * 4 or 1 

ax 


3T (32) 

3T k z 

s * z 2 + (y - >v) 2 ’ 

1^2 = n-j k * 3 or 2 

3T< 32) -(y - \) 

31 z 2 + (y -\) 2 

(2.6.23) 

— — * unchanged 

- 

jtJ 41 * 

— — = unchanged 

k = 4 or 1 

n' 4 " 

— * — * unchanged 
aZ 


All of the quantities of Appendix D are 
values, except that Hq 2 is replaced by 

now recalculated using these modified 

II m ** ’ - * 

‘* 4lY ’ ,(41) *” (2 - 6 - 241 

h o2 ”»i 77^ 

it + 

The induced velocities from the last wake element are added to the other dipole 
velocities of the lifting strip in the ordinary way. 


2.6.4 Some Special Situations 

Two special situations exist where panels must be placed inside the body 
surf.ce. No normal-velocity boundary condition can be applied at such elements 
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and no source density should be applied to them. However, they do have vor- 
ticity and this must be accounted for. 

The first situation occurs when a portion of the bod, intersects a lifting 
portion at a finite angle (often nearly nonnal) without breaking the continui 
If the trailing edge. An exauple is provided by the wing-pylon '" tersect 
shown in Fig. 5. A certain portion of the lifting body surface is inside 
pylon. However, the underlying dipole distribution should be continuous 
through this region to avoid nu-erical difficulties. Thus, asfar as vorticity 
calculations are concerned, the ‘inside- panels are normal ueubers of the 
ing strips to which they belong. But they are Ignored as far as source calcu- 
lations or boundary conditions are concerned. Such panels are designated 
"Ignored panels.- They usually co^wise only part of a lifting strip. 

The second situation occurs when a lifting portion of the 
another portion at a finite angle (often nearly normal). The inportant case 
of this is the wing-fuselage Intersection, as illustrated in Fig. • * ' 

well-known, the local -section lift coefficient- on the wing does not fall to 
zero at the fuselage intersection. Thus, the underlying dipole strength on the 
11-line lying along the intersection is not zero. However, t e ' 
cannot singly be terminated, because that would result 1. . . 
vortex filament right on the surface. Accordingly, an additional or extra 
lifting strip is added to the lifting section (see Fig. 6). It is either 
Hrst or Tlast strip of the lifting section. The extra strip lies inside 
the other body and is a complete lifting strip including wake. No source A n- 
sities or normal-velocity boundary conditions are applied to the panels of 
extra strip. The underlying dipole strength is taken constant in the span- 
wise- direction across the extra strip. The value of the dipole strength on 
the extra strip has nonzero dipole strength and may lead to a concentrated edge 
vortex in the streamwise direction. For example, as shown in Fig. b, 
vortex may lie along the fuselage centerline and its downstream extensor., 
the lifting configuration has a right-and-left symmetry, e.g.. a fu«lagew,th 
noth wings, and if the flow is also symmetric, e.g. zero yaw. the extra P 
for the right and left sides have the same strengths on their interior e ges. 
Thus, in this case the edge vortices cancel. If, however, the lift is not 
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symmetric, there will be an edge vortex. This is unavoidable because it is 
physically real. An example is the hub vortex of a propeller. This also 
occurs at a tip tank, which is essentially a small fuselage with only one wing. 

2.6.5 Assembly of the Vorticity Onset Flows 

As described in Sections (2.6.1) and (2.6.2), the velocity induced by the vor- 
ticity on a panel and the associated edge vortices fall naturally into two 
parts - one proportional to the value of B on the first N-line and one propor- 
tional to the value of B on the second N-line. These are summed over the 
lifting strip to yield two vorticity onset flows for each lifting strip. In 
general, each onset flow has three components at every control point. Specif- 
icity, 

♦(F) . str JP k #<F> 

T ik ' J ¥ ij 

k * 1, 2 L (2.6.25) 

W - st ? k *S» 

where L is the number of lifting strips. The summations of Eq. (2.6.25) are 
over a complete lifting strip including the wake elements. If a lifting sec- 
tion begins with an "extra strip" (Section 2.6.4), both velocities and fjjp 
for the extra strip are added to the velocity corresponding to the first 
ordinary strip of the section. Similarly, if the last strip of a lifting 
section is an extra strip, both velocities for the extra strip are added to 
the of the last ordinary lifting strip of the section. (This gives an 
underlying dipole strength on the extra strip that is constant at a value equal 
to that attained on the adjacent lifting strip along the camm>n N-line of vhe 
two strips.) Thus, the calculation of Eq. (2.6.25) gives an N x L matrix of 
velocities at the control points, where L refers to ordinary lifting strips 
only. Since L is small conpared to N, these matrices are small compared to the 
source-velocity matrices. Each of the velocities, Eq. (2.6.25), represents the 
velocity due to an underlying dipole distribution of the strip that has slope 
unity on one N-line and zero on the other with a linear "spanwise" variation 
in between. 
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The characteristic onset flow velocities due to a strip are 


C * 4 s 1 * 4' 

1 14’ - 


(2.6.26) 


The first velocity of Eq. (2.6.26) is that due to an underlying dipole distri- 
bution on the strip that is constant in the "spanwise“ direction. The second 
velocity is that due to a dipole distribution that varies linearly in the 
"spanwise" direction and has zero value at "midspan." These velocities are 
used to form the basic circulatory onset flows 


If the "step function" option for bound vorticlty Is used, the proper form of 
the dipole distribution is simply constant in the "spanwise" direction over 
a lifting strip, and the velocity is precisely the onset flow. Thus, for 
this option, the vorticity onset flows are 

*< k) = k = 1, 2, .... L (2.6.27) 

The above yields L onset flows, each of which corresponds to a unit value of 
the "streawwise" dipole derivative B on one lifting strip and zero values of B 
on all other lifting strips. 


The machinery for the "piecewise linear" option for bound vorticity is somewhat 
more complicated. The "spanwise" variation of the "streanwise" dipole deriva- 
tive B (bound vorticity) is linear in the "spanwise" direction across a lifting 
strip. Thus, the velocity at the i-th point (control point or off-body point) 
due to the k-th strip is 

f, (“HP k) = 5'°' D k + ( 2 . 6 . 28 ) 

where w^ is the "spanwise" width of the strip. B' is the "spanwise" deriva- 
tive of B, and subscripts k denote quantities associated with the k-th lifting 
strip. The derivative V k is evaluated by a parabolic fit through B J( _ 1 , 

B k and B k+ j. Specifically, define 
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\ * Vi 


D k ’ • W k 7 T7JTV 


+ W, 


k+1 


"k + "k-1 


E k * » k + T77T5Q - - k+ ,l 


w.. 


[ 


w k * w k+1 
\ + W k-1 


B "k r VVl ] 

F k *k + w k+l 

Then Eq. (2.6.28) is approxi*ated numerically by 


*k + w k-l 

\ + Vl 


(2.6.29) 


f, (strip k) • I»\ * 1 \ 1 ] * \\ * Wl 1 


(2.6.30) 


The velocity 6). (2.6.30) contains v.loes of the -st™»1se- dipole dwlvetlve 
B for three consecutive strips. However, a proper circulation onset flou Is 
proportional to the value of * on a single strip. Since each Bk «*•” 

5. ( s tr1p k) for three consecutive strips. Its three contributions My he 
sunned to give the basic vorticity onset flow. 

f< k) ■ C ♦ t * ^ ,£ k ♦ C.v, l2 • 6 • 31, 

in performing the above parabolic fit Eg. (2.6.30). the values of tbe functlon 
B to be fit are. of course, the values of bound vorticity on the strips. Uc 
of these has been associated with an abscissa or independent .arlable that 
expresses the spanwlse position of each strip. Differences of these abscissas 
appear as cognations of the widths v Calculation of the w^ is "»t 
„b!tous. because In genera, the -span- or width of each strip is no constant 
but varies in the -chordwlse’ direction. An average width is calculated for 
each strip and used in the calculation above. 

The calculation.! machinery of the pro*., insures that the underlying dipole 
strength varies linearly In arc length along an N-llne, l.e. that dipole 
strength equals Bl, where B Is a constant to be determined .ne t Is arc 
length measured fro. the lower surface trailing edge. In particular, at the 
upper surface trailing edge, the dipole strength Is BL (total), flils is th 
circulation about the B-llne and the value of vorticity that carries into the 
wake. The machinery above fits the spwwise distribution of dipole derivative 


\ 


B, but It nukes better sense physically to fit the circulation distribution BL 
(total). This is a smoother function because it is independent of planform 

breaks. 

The code has the option of fitting either B or BL (total). While the above 
concept is somewhat complicated to explain, its numerical implementation is 
simplicity itself. All that is necessary is to divide the vorticity onset 
flows associated with each N-line by the corresponding values of L (total). 

Thei. the values of "B" that are solved for will really be values of BL (total). 
Thus in Eq. (2.6.26) 

a(F) 

LJtotaTT replaces V ik 

F (2.6.32) 

qfeuTT r *>”* ce5 *» 

No other changes are necessary. These are added to produce a single dipole 
onset flow per strip and a complete flow solution obtained for it. 


2.7 The Kutta Condition 


The Kutta condition is applied as a condition of pressure equality on the upper 
and lower surfaces of the trailing edge, which amounts to a condition of equal 
velocity magnitude. As a numerical approximation, the Kutta condition may be 
applied by equating pressures at the control points of the two panels adjacent 
to the trailing edge on the upper and lower surfaces of the wing. Alterna- 
tively, velocities at the upper surface control points of the few panels near- 
est the trailing edge can be extrapolated to obtain velocity components "at" 
the trailing-edge upper surface, and the same could be done for the lower 
surface. This last allows application of the Kutta condition more nearly at 
the trailing edge, and the analogy of this procedure yields considerable 
ingjrovement in accuracy in two-dimensional cases. This is the option used in 

the present method. 
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However the point of application of the Kutta condition is chosen, the logic 
of the calculation is the same. In particular, a velocity vector can be calc- 
ulated at the upper and at the lower trail ing-edge point for each strip of 
panels (Fig. 1). From the discussion of Sections 2.5 and 2.6.5 it is clear 
that the velocity vector at the i-th control point is given by 



N „ 

I V^o 
j=l 


L 

i i + I 
ij J k=l 



(2.7.1) 


where L is the onset flow (usually a uniform stream). Initially, the 7j and 
B. are unknown, but it can be seen that depends on them linearly. If 
velocities are extrapolated to the trailing edge, this linear dependence 
remains. Let superscripts U and L denote velocities at the upper and lower 
trailing edge, respectively. Further let subscript m - 1, 2, .... L denote 

conditions on the m-th lifting strip of panels. Thus vi U) 4 U 

the velocity vectors at the upper and lower trailing edge, respectively, of the 

„_th lifting strip. The condition that these two velocity vectors have equal 

magnitu des may be written In terms of dot products as 

fill . Iff' ♦ r U J . » ,U • lt‘ w * tf-’j (2.7.2) 

Apply)*, Eq. (2.7.2) .t tlw tr.111.9 «He Of e«A Hftlnq strip yields t qued- 
ratic equations in the (N + L) unknowns Oj and 


2.8 Iterative Matrix Solution 


The velocity induced at the 1-th control point by the source and vortex 
singularities is given in Eq. (2.7.1). Taking the scalar products of these 
velocity vectors with the unit normal vector ^ of each control point gives 
an N x N matrix of source influence coefficients and an N x L matrix of vortex 
influence coefficients defining the normal-velocity Influence of each element- 
ary singularity distribution at every control point. If we define the solution 
vector X whose entries are the source strengths U = 1. H) followed by 
the vortex strengths B*. (k » 1, L>, the condition of zero normal velocity 

on the body can be written in the form 
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AX = R 


( 2 . 8 . 1 ) 


where A is the N x (N+L) matrix of source and vortex norml velocity influence 
coefficients, and the right side in Eq. (2.8.1) is the negative of the normal 
component of the freestream velocity. It is important to note that A is pure 
geometric and does not depend on the onset flow, which enters only the r,ght 

side. 

Equation (2.8.1) defines a system of H linear equations in the (Net) unknown 
singularity strengths. A further set of l equations therefore 
complete the fonnulation of the problem. These equations are provided by the 
Kutta condition, whose derivation is outlined in the previous section. This 
condition ensures that the co^uted upper and lower surface P r «sures mac 
the trailing edge. The resulting equation, (2.7.2), can be written in the for. 

,»<«> - f (L >) . f a 


( 2 . 8 . 2 ) 


av 


where I* 1 * and f (U are the upper and lower surface trailing edge velocities 
while f is the average of these two velocities. 

One Kutt. condition is applied for each lifting strip, giving a set of » linear 
and L nonlinear .quations to be solved for the H+L unknowns. For caplex con- 
figurations, H can be large, up to 2000 in the current version of the code, 
while L is typically between one and two orders of magnitude less. 

The costing tie* required for the solution of the linear equations by a 
direct solution is proportional to while that required for an iterative 
matrix solution is proportional to N 2 per iteration. Therefore provide, 
that the nuafcer of iterations required to obtain a converged solutions is rel- 
atively small (cohered with the n«*er of unknowns), there ,s a large benef. 
to be obtained through the use of an iterative matrix solution. The sdieme 
adopted here is an accelerated block Gauss-Siedel iterative procedure which has 
been shown to give rapidly convergent solutions for a wide rag. of g~«tries 
(p e f 8 ) This section will outline the details of this Iterative scheme, ad 
Append^x H will present the details of the acceleration schae which has bea 


3495H 


27 


adopted in order to improve the speed and the stability of the convergence 
procedure . 


As pointed out above, the Kutta condition to be applied is nonlinear, and so 
it must be linearized in some manner consistent with the iterative solution 
procedure which is to be applied. If we introduce the subscript K to denote 
quantities evaluated after the K-th iteration, Eq. (2.8.2) can be written as 

(» <U) - » (U ) |C . - 0 (2.8.3) 

so that the average velocity from the previous iteration is used. For the 
first iteration, $ ay is replaced by a vector along the local trai ling-edge 
bisector, which ensures that the physically meaningful root to the Kutta con- 
dition is selected in which both the upper and lower surface velocities leave 
the body at the trailing edge. 

2.8.1 Block Gauss-Si edel Iterative Scheme 

It was shown by Hess (Ref. 9) that the Gauss-Siedel iterative matrix solution 
scheme converges very rapidly for simple external flow problems. This 
approach, which is described by Varga (Ref. 10), relies on solving a succes- 
sion of lower triangular matrix equations of the fora 

AY = R - A X _ (2.8.4) 

A JTK K *u A K-1 

where X K is the K-th approximation to the solution. The matrix represents 
the diagonal and lower triangular part of A, while A y represents the upper 
triangular part. 

For lifting flow problems there «s a strong coupling between the source and 
dipole strengths for a given lifting strip. Therefore, in order to maintain 
the diagonal dominance of the matrix, it is necessary to adopt a block Gauss- 
Siedel scheme. The particular approach used here takes the source strengths and 
the associated dipole strength for each lifting strip as separate blocks in the 
solution vector. In this way the normal velocity conditions for a given lift- 
ing strip are satisfied simultaneously along with the Kutta condition before 
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proceeding with the solution for the next block. This is not the only way in 
which the block structure could be implemented. Reference 11, for instance, 
groups all of the dipole unknowns together as a single block in the matrix. 
However, for the nonlinear Kutta condition, the approach adopted here is more 
convenient. For nonlifting sections of a configuration, the choice of the 
block structure is less crucial. As the block size for such panels is 
increased, the computational effort is increased but the rate of convergence 
is also increased. The 'ise of a block size of 50 has been found to give a good 
compromise. 

This iterative scheme is equivalent to solving a series of quasi-two- 
dimensional problems corresponding to each block in the matrix. The onset 
flow for each of these calculations includes the current effects of all the 
other panels on the body. Therefore, as the solution converges, a fully 
consistent three-dimensional solution is obtained. 

The iterative solution procedure can be broken down into two steps, the first 
of which involves the calculation of the rl^it-hand side of Eg. (2.8.4) based 
on the previous solution. 


RHS K _1 = R " A u X k1 (2.8.5) 

The second step is the calculation of the new approximation 

= RHS^_j (2.8.6) 

Each of these steps is performed successively for each block of unknown source 
strengths, each of which involves the direct solution of a small set of simul- 
taneous equations. In addition, for lifting strips, the dipole strength is 
computed by satisfying the Kutta condition for that particular strip. 

At this stage, »t should be pointed out that, for large matrix equations, the 
whole coefficient matrix A cannot be stored in the computer memory at one time. 
The order in which the matrix is formed and stored on disc will therefore 
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• +h» utrix solution scheme is formulated. In gen- 

influence the way in which the matrix s matrix 

ite rr c r * ». »*.. 

ZVlZt schemecan be to matrices store, either t, r«s or h, 

col-s. The Matrix multiplication operations Involve. » bo « - 

xno (2 8 6) can be acco.pl i shed for either row or column stored matrices. 

arises from the order in which the multiplication loops are 

nested. 

2.8.2 Convergence Acceleration Scheme 

.. ^ ir rsm rrr;;« 

... rr: “r 

r«mr rjrss. 

U^ift systems or win^pylon/n«*iT« -runt*.. «• 

^ a . rases it can fail to converge entirely. 

scheme becomes worse, tedutioues are based on the extrapolation 

Traditional convergence acceleration techniques are oasea 

of the solution m*1ng use of the asymptotic converge. "“•***£“ 

, discussion of several iterative schemes. recommending compo 

relaxation to d-P out any -cl^ons^n 

for cm ”': ::zTrZ • 

unsuitable. new — 

Much can be app le a sdl _. h as been found to give improved converg- 

r^rrrs^ .is. -»•*. " c ^ to 0 f be 

obtained for cases which are well outside the normal ran^ of c.nv«-genc 
the basic Gauss-Siedel iterative scheme. 
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The scheme adopted is a relaxation method in which, after each iteration, an 
improved solution is defined as a linear animation of the earlier approxima- 
tions. However, the relaxation coefficients are computed after each iteration 
in such a way that the residual error of the new approximation is minimized. 

To implement the scheme described in the previous paragraph, it is necessary 
to define the residual vector after each iteration. For the linear normal 
velocity equations, this residual is defined by 

RES k = R - AX k (2.8.7) 

It is clearly undesirable to have to evaluate this expression after each iter- 
ation since this would involve (N 2 ) operations which is equivalent to an 
additional Iteration. However, by separating the matrix A into its triangular 
parts, and applying Eqs. (2.8.5) and (2.8.6), the residual vector is given by 

res k = rhs k - rhs k . 1 (2 ‘ 8 * 8) 

and this expression can be easily evaluated at the end of each Iteration. The 
residual for the Kutta conditions can be evaluated by computing the trail mg- 
edge velocities after each iteration. Substitution in Eq. (2.8.3) then gives a 
value for the Kutta residual for each unknown dipole strength. 

Apart from the calculation of the residual vector, defined by Eq. (2.8.7), the 
convergence acceleration scheme presented here does not depend in any way on 
the details of the Gauss-Siedel iteration. The schemes are applied as two 
independent steps of the overall iterative procedure. In the following outline 
we will define RES K to be the residual vector, including both the normal 
velocity and the Kutta condition residuals after the K-th iteration correspond- 
ing to the solution vector X K , which also includes both the source and the 

dipole unknowns. 

Given a set of approximations to the solution, X Q , Xj, ...» X K » we can 
define a new approximation by 
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(2.8.9) 


K 

X' » l Xifi 
i=0 

where f Q , fj, ..., f k are thr acceleration coefficients which are yet to 
be determined. In general, the first approximation, X Q , which is the start- 
ing solution to the iterative procedure, is taken to be the zero vector. 
Therefore, Eq. (2.8.9) defines a set of K independent approximations to the 
solution vector, and so it is convenient to constrain the acceleration 
coefficients so that 


f fi = 1 (2.8.10) 

i=0 

It now follows from Eq. (2.8.7) that the new residual vector is given by 

K 

RES* * l RESffi (2.8.11) 

1*0 

It should be noted that, while this equation Is exact for the normal velocity 
residuals which satisfy a linear equation, it is only approximate for the 
nonlinear Kutta residuals. However, this approximation is consistent with the 
linearization applied in the solution of the Kutta condition, and it is a good 
approximation for this application. 

As the coefficients f i vary, Eq. (2.8.9) defines a family of approximations 
to the solution of both Eqs. (2.8.1) and (2.8.2) for which the corresponding 
residual is given by Eq. (2.8.11). In order to minimize the error for this new 
solution, a single scalar measure of the error is required. The sum of the 
squares of the opponents of the residual vector provides a suitable error 
measure. In matrix notation this quantity can be evaluated in terms of the 
non of the residual vector which is defined by 

| IRES' | I 2 » [RES'] T RES' (2.8.12) 

where [RES'] 1 is the matrix transpose of RES'. 
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This equation defines a quadratic function of each of the variables f.. 

Since this function is non-negative, it follows that its minimum value must 
occur at the point at which 

Si (RES' | | 2 /9f i =0 for i = 0, 1, .... K-l (2.8.13) 

This provides a set of K linear equations which, together with Eq. (2.8.10), 
can be used to determine the acceleration coefficients completely. 

Full details of the derivation of this set of equations and their solution are 
given in Appendix M. 

This acceleration scheme involves two principal computational tasks. The first 
is the calculation of the acceleration coefficients, which in turn involves the 
calculation of the scalar products between every pair of residual vectors. The 
second is the application of these coefficients to the calculation of an 
improved solution and its corresponding residual vector. Both of these tasks 
will involve on the order of (KM) operations while each iteration requires N 2 
operations. Therefore, since K, the ntafcer of previous solutions, i: very 
much less than N, the additional computation Introduced by the acceleration 
scheme is small. However, it has been found to have a significant effect on 
the rate of convergence of the scheme, while also giving an improved stability 
enabling converged solutions to be obtained for cases which are well outside 
the normal range of convergence of the basic Gauss-Siedel iterative scheme. 
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3.0 THE I MET PROCEDURE 


3.1 General Description 


The inlet procedure employs the above-described panel method to calculate 

ob^the f l' 0 " “ ^ tlW in ' et “ h,Ch ^ the " Umr '> “«*><"«< to 
may be 1! ! < ‘ eSl ’' ed CI, " d1tfo "- Specifically, solutions 

may be obtained for any angle of attack or yaw, Mach number, and mass flow 

e. The computat.onal effort required to perform the combination for a par- 
a ar operating condition is a small fraction of that required for the i„U- 

of olrat ! the fmd31KnU ' S0,utiMS - solutions for any number 

t rr: 9 T my be as needed, at any time 

after the fundamental solutions have been calculated. 

The nwerical efficiency of this inlet procedure is realized because the funda- 

r™ ZT mS " 0kU,ned f0r now. and then c^ii T 

«n^tod for compressibility effects. * key «l««nt In this apprsuch Is an 
accurate and general cmpressibility correction that may be applied to the 
«c-press1bl. n«. about the sa. Inlet, as opposed to the stands Soethert 
P r * ,l,re * * ^ooh-romdwr -dependent stretched version of the 

The compressibility correction u*4 1, the Liebl.ln and Stock«n 
method. Ref. 13, which is described In Appendix K. This procedure has been 

n«,T d « y C T r, “" * Uh ' Xperi,,enU1 “***' «<rfs. 13-15. For Internal 

oener 2 ! **" SU| ’* rS0, ’ 1c f, « “^"out shocks, and it has been 

generalized to external flow about wings. Ref. 76. 

fromthebeginningiRef. „ this work has had t« principal aims: co*uta- 
t.onal efficiency for arbitrary geometries, which is discussed above, and user 
Mentation, which has been obtained principally by including a nutter of 
graphical output features. The main capabilities of the programs of Refs. 1 
an are surface streamline tracing and isoplotting of various flow quantities 
both on the surface and over cross sections. Both of these have been inroved 
by providing the caprtlHty of drawing curves across section boundaries. That 
IS, the panels may be grouped into logically independent but physically contig- 
uous sections or networks, and the plotting routine can draw streamlines or 
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i socurves across these boundaries. A major new graphical feature is the por- 
trayal of the surface or off-body velocity field by means of a set of vectors 
having the velocity magnitude and direction at all points. This type of 
picture has proven very useful in applications. 

3.2 The Fundamental Flow Solutions 

First, the definition of a flow solution must be described. In the present 
context these are incompressible flows. Every flow solution corresponds to a 
certain -onset flow- which is the flow incident to the body. In general the 
onset flow satisfies neither the normal-velocity boundary conditions nor the 
Kutta conditions. The source densities o^ and the dipole derivatives By 
(bound vorticity strengths) are adjusted to satisfy these condltions The most 
comnon onset flow is a uniform stream, but as will be seen, other onset f ows 
are also necessary. For this reason, the onset flow vector at the panel con- 
trol points is written » 0 , to slum that it may vary fro. point to point. Then 

the velocity at the 1-th control point Is 


V i 


H * 

l *1j a j + 

j=l 1J J 


I ^By 

k=l 1 K 


oi 


( 3 . 2 . 1 ) 


This replaces Eg. (2.7.1), and the method of Section 2.7 and Section 2.8 give 
the values of o. and B. corresponding to that particular onset flow. 

When these values are Inserted into Eg. (3.2.1) and the indicated sumptions 
performed, the resulting ->t of is designated a flow solution. 

The set of fundamental flow solutions that are superposed by the combination 
program to obtain flow about the inlet at arbitrarily prescribed operating 
conditions may be described most easily in terms of two types of flow. The 
first is flow about the inlet due to a unit frees .ream at prescribed ang.e o 
attack and yaw with no effort to control mass flow through the inlet. In the 
nonlifting nethods of Refs. 1 and 2, there were always three such flows: zero 

angle of attack and yaw, 90- angle of attack, and 90- angle of yaw. In lifting 
cases, however, the latter two flows make no sense. If the inlet has a 
leading-edge slat, for some circuuferential locations, the trailing edge is 
the upstream point of the airfoil section. The result can be nonconvergence 
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of the iterative solution. In the present program a number of angle-of -attack/ 
yaw combinations are input by the user. It is preferable to choose these com- 
binations in the range where the user's interest ultimately will lie. The 
other type of solution is the static, for which the inlet ingests fluid and the 
flow is quiescent at infinity. The inlet methods of Refs. 1 and 2 use differ- 
ent mathematical devices to produce the static solution. Reference 1 uses a 
constant vorticity distribution over the inlet surface, as illustrated in Fig. 
7a. This has some features in common with the surface vorticity used on the 
slats to generate lift, but is also has several differences. No Kutta condi- 
tion is applied on the inlet, and no distributio'i of vorticity is solved for. 
Instead a single parameter, total vorticity strength, is adjusted to satisfy a 
single condition, mass flow through the inlet. If auxiliary inlets are pres- 
ent, the topology of the configuration does not permit use of surface vortic- 
ity. Accordingly, In the method of Ref. 2, the mechanism of the static solu- 
tion is a single ring vortex located well downstream in the Inlet, as shown in 
Fig. 7b. The strip vorticity option of Fig. 7a gives a superior static solu- 
tion, and It is used In all cases except for the Infrequently occurring one 
where an auxiliary inlet Is present. In the very Infrequent case where there 
are two independent mass flow rates, e.g. an "Inlet within an inlet," the above 
mechanisms have to be applied to each Inlet separately. 

Because of the wide variety of cases to which the present method may be 
applied, some flexibility is necessary in the choice of fundamental flow solu- 
tions. For example, while the static solution has a sensible Kutta condition 
for an inlet with leading-edge slat, the same probably cannot be said for an 
inlet on a wing, where the static flow near the wing trailing edge is more-or- 
less parallel to it. Similarly, the high Inclination angles at which a slatted 
inlet can operate at high mass flow rates lead to the above-described diffi- 
culty for 90° if no mass flow control Is exercised. Thus, in general the 
fundamental flow solutions should all contain combinations of an inclined 
freestream and a static condition. This is perfectly permissible as long as 
the flow solutions contain all the independent possibilities, e.g. at least 
two angles of attack, yaw, and mass flow rate. 
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3.3 The Combination Program 

The fundamental flow solutions and the body geometry are accessed by the com- 
bination program. At this stage also are input any off-body points and inlet 
cross-sections where the flow output is desired. A cross section is a panel 
network extending across the interior of the inlet. Flow quantities are com- 
puted at panel centers and a total mass flux for the cross section as a whole 
is evaluated. One cross section is designated the control station, and it is 
there that the mass flow condition is applied. In preparation for this, the 
average velocity at the control station V is computed for each fundamental 
solution. 

The flow condition input to the combination program consists of flow conditions 
at infinity and at the control station. The various possibilities are pre- 
sented in Appendix 0. The key quantity in the ctmbination is the equivalent 
incompressible velocity, which is denoted with a prime. In particular, is 
the equivalent incompressible freestream velocity (Eq. CO-13)} and ¥ c is the 
equivalent incompressible average velocity at the control station (Eq. (0.20)). 
In all cases V equals Y multiplied by the local static-to-total density ratio, 
and the flow direction is unchanged. 

In order to compute the combined flow for a given set of flow conditions, a 
number of the fundamental flows are combined linearly. In general, three lin- 
early independent fundamental flows are required to satisfy the conditions at 
infinity while an additional static solution solution is required for each 
independent mass-flow condition. However, for flows without yaw, the number 
of fundamental flows required is reduced by one. In the fundamental solution 
mode, a nunber of user-specified fundamental solutions are obtained including 
at least one yaw solution if combined solutions with yaw will be required. The 
range of angles of attack and yaw specified should preferably span the complete 
range of confined solutions which will be required. When the combination pro- 
gram is run, the code will automatically select the closest linearly independ- 
ent solutions to be used for the combination. This procedure is required by 
the nonlinearity in the potential flow solution which is introduced by the 
Kutta condition. 
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To Illustrate the combination procedure, consider the case In which four 
fundamental solutions have been selected by the code. Let these Individual 
solutions be denoted by the superscript m, and let a^ represent the unknown 
combination constants. The equivalent Incompressible velocity for the combined 
flow is 

p. = l a (3.3.1) 

1 ra=l 1 

where the combination constants are Initially unkmwn. Meeting the pre- 
scribed flow conditions at infinity and at the control station requires 



This defines four equations (one vector, one scalar) for the four unknown a R . 
Once computed they are inserted In Eq. (3.3.1) to obtain t^ f tdtlch Is used 
In the compressibility correction (Appendix A). 


4.0 CALCULATED RESULTS 


The Method described here win be Illustrated using three separate geometrical 
configurations. The first two cases represent complex three-dimensional flows 
Involving inlets with lifting slats on the leading edge of the cowl. The third 
configuration is a simple nonlifting axi symmetric inlet which illustrates some 
improvement in the computed results, as compared with the nonlifting method 
presented in Ref. 2. 

The first geometry discussed Is an axisymetrlc inlet with a centerbody and 
leading-edge slat shown in Fig. 9. The cross-section shown In Fig. 9b illus- 
trates the relationship of the leading-edge slat to the cowl. The ability of 
colored shaded graphics to portray complex three-dimensional bodies is illus- 
trated in Figs. 9c and 9d. Figure 10 shows a comparison of the current method 
with the axl symmetric method of Ref. 18 for a cotfrined flow along the axis of 
the Inlet with an average velocity on the fan face, x - Z.S3, of twice the 
frees tream. A surface vortlclty distribution on the cowl was used to generate 
the static solution in both the axl symmetric and the three-dimensional calcu- 
lations. The pressure distributions on the cowl and the centerbody, shown In 
Fig. 10a, agree very closely while Fig. 10b shows that on the slat there is a 
small difference In the leading-edge pressure peak when compared with the 
axlsjmuetrlc result. 

The remaining results for this configuration are presented for three incompres- 
sible flow conditions. The first is a pure static flow with no flow at infin- 
ity, and the second is a -pure freestream- solution at zero angles of attack 
and yaw, with no surface vortlclty on the cowl. The third solution Is a com- 
bined flow at 40° angle of attack, zero yaw and an average fan face velocity 
twice the freestream velocity. Figures 11-13 illustrate the flowfleld across 
the inlet for these three flow conditions, the vectors drawn being proportional 
to the local flow velocity. The flowfield velocity vectors are computed in the 
plane through the inlet axis tilted 15° from the center plane of the inlet. 

The boundary lines shown on these figures represent the boundaries of the off- 
body flowfleld rakes used to compute the flowfleld rather than the exact 
aerodynamic surfaces. For clarity, the approximate body locations are shaded. 
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The static case shown In Fig* " 111 *tr.t.s the rapid <*xease In reloclt, 
magnitude ahead of the Inlet, uhlle the expanded view of ** *'* ^’"\ 
reveals a local flow environment similar to that of a conventional flap with 
the upstream stagnation point occurring In the vicinity of the slat leading 
edge. The axi symmetric freestream condition shown in Fig. gives 
unconventional flowfleld. In this case the -upstream- attachment line on t 
slat occurs on the forward-facing surface of the slat closer t 
edge than the leading edge, while the attachment line on the cowl occurs vir- 
tually under the leading edge of the slat. Figure 13 ’''“ttrates » 
now in which the freestream is at 40- to the inlet axis while the -nter a 
velocity is twice the freest ream. In this case the flow .*e d 

on the centerbody Is very different than that 

the vicinity of the slat, the flowfleld Is qualitatively the same. 


Figures 14-H pro** the canted flowfleld 1t*» 

conditions 1m the smme Off -body plane. Figure 14 

nat— f the slat ftamfleld 1. the stetlc case with the «p1d 

.. f v^ neimr nmni tilt Inner lip of tilt cowl# 0t tiw 
tlon occurring ts the flow gets mma me inner r . 

other hand. Figs. 15 and 16 show the extreme pressure gradients occurring on 
the slat where the flow twms around the leading edge. 

The second configuration considered Is essentially the previous configuration 
but further co, Heated by the addition of another loading-edge slat, this ^ 
only a part-di ff ere n ce slat, however. » section t rough ... 

the configuration is shown in Fig. 17. Results are presen or 
uration at a fblned solution of zero yaw, 40- angle of attack and 
velocity of twice freestrean. 

Figure 18 compares the co*> uted pressure distribution, plotted against radial 
distance, on the -in slat at three different 

the corresponding results counted for the single ^ 

top of the inlet the presence of the auxIIW *1* **" ' "Hot 

effect on the pressure. However, 1 u the z > 0 plane there Is gn ^ 

reduction in the leading-edge presswe peak, —lie c ose 0 

inlet the second slat greatly reduces both the leading and trailing edge 
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pressure peaks on the main slat. Figure 19 shows the surface velocity distri- 
bution on both the Inner and outer cowl surfaces for the same flow conditions. 
It can be seen fro* Fig. 19b that the nonaxlal nature of the flow persists 
throughout even the interior of the inlet. This is presumably due to the 
presence of the ingested tip vortices trailing from each end of the part- 
circumferential slat which will induce some swirl into the internal flow. 

The third configuration considered is a single 72-panel (on th* half-body") 
round inlet, as shown in Fig. 20. This simple geometry was used to demonstrate 
the improvement gained in the new source derivative fitting algorithm in use 
in the present code. Figures 2la and b present the variation in peak Cp versus 
theta (measured circumferentially), and the variation in Cp versus axial dis- 
tance at a fixed theta value (0 * 75°), respectively. The results are com- 
pared with those obtained by the method described in Ref. 2 in which a "least- 
squares" fitting procedure was used to compute the source derivative effects. 
The present approach, described in Section 2.5.1, demonstrates the Improved 
level of ax i symmetry which is obtained by the new formulation. 
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5.0 INPUT INSTRUCTIONS FOR THE HIGHER-ORDER POTENTIAL FLCW PROGRAM (DF12) 

5.1 Introduction to the System 

The computer code is actually a collection of pre- and post -processing programs 
grouped around the potential-flow program. It can be thought of as a system 
of programs designed to "talk" to each other via saved datasets. These pro- 
grams are: 

1 . PRE-PROCESSOR: parametric cubics patch fitting of 3-D coordinate data. 

2. FUNDAMENTAL POTENTIAL FLOW SOLVER: (DF12: Mode 1). 

3. COMBINATION OF FUNDAMENTAL FLOWS (DF12: Mode 2). 

4. POST-PROCESSOR: ISOPLOT - plots iso-contours (on- or off-body). 

5. POST-PROCESSOR: VECPLOT - plots velocity vectors (on- or off-body). 

6. POST-PROCESSOR: TRACE-ON - calculates streamlines (on-body only). 


Operation of these codes is facilitated by a set of Interactive "submit CLISTS" 
and associated FORTRAN programs. A single CLIST controls the operation of 
programs 1, 2 and 3; separate CLISTS exist for each of programs 4, 5 and 6. 
While all these CLISTS have been designed for an IBM mainframe running TSO In 
an MYS/XA environment, similar interactive submittal procedures can easily be 
written for other systems to accomplish the same purpose, viz. simplify the 
user's job of running cases and enhance his ability to investigate both the 
quality and significance of the computed results. 

5.2 Discussion of the Individual Programs 

Since several of the programs can "communicate" with each other via saved data- 
sets, a great deal of flexibility exists concerning the sequence in which the 
programs may be executed. For example. Program 1 can talk to either Programs 
2 or 3, but is not always required to run Programs 2 or 3. Programs 4, 5 and 6 
can talk to Programs 2 or 3, but only if the appropriate dataset from 2 or 3 
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was saved. To understand the possible Interactions between programs. It Is 
best to consider each one separately, first. 

5.2.1 The PRE-PROCESSOR: PC-PATCH 


This program is designed to take a user-defined set of 3-D Cartesian coord- 
inates and fit a set of parametric bi-cubic patches to the implied surface. 
The input consists of a formatted "card-image" (i.e. 80-column, fixed block) 
dataset which contains the corner points of panels distributed on the surface 
of the body about which the flow is to be calculated (see Appendix P)„ The 
format of this data is: 


cc 1-10 X 
cc 11-20 Y 
cc 21-30 Z 

cc 31 ISTAT 


cc 32 LABEL 


cc 33 MCURY 


Cartesian coordinates 


0 3 this point Is on the same (Mine as previous point 

1 3 this point starts a new N-llne 

2 ■ this point starts a new section 

0 - this is an NLIF (nonlifting) section 

1 3 this Is a LIFT (lifting) section (I.e. has Kutta 

condition) 

(Note: all LIFT sections. If any, must precede all 

other section types.) 

2 3 this is a HAKE section (all HAKE sections. If any, 

must come last on the Input geometry dataset, after 

all other OREL =0,1, 3 and 4 sections) 

3 3 this is a DBLT (doublet) section 

4 3 this is a SRFV (surface vorticity) section 

(Note: may not use both DBLT and SRFV sections at 
the same time In a Node 1 case.) 

5 3 this is a FLUX section (allowed as input to Node 2 

cases only) 

0 = automatic M-line curvature selection (curved unless 

LABEL 3 ! ) 

1 = H-lines are all curved (even if LABEL 3 !) 


Note that LABEL and MCURV only apply to ISTAT=2 points, i.e. they need be 
entered only once on a section. (For a discussion of the limits on the numbers 
of points, panels, etc., see Section 6.2.4, Program Limits.) 


These panel coordinate data are "fit" with parametric cubic patches (see 
Appendix B) and written to an output dataset hereafter referred to as a "PCU"- 
dataset (Parametric Cubics Unformatted). The PCU dataset serves as a true 
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surface definition for the higher-order potential -flow solver, but It Is not 
required that the pre-processor supplied with this system be used to generate 
that PCU-dataset: any "PC-fitting prograa may be used, as long as the follow 


ing PCU-"format" 

is observed: 


Record #1: 

I FORM 

A single integer (use "1") specifying the PCU- 
format 

Record #2: 

NSECT 

NPATT 

NTYPE(6) 

HEAD(9) 

Number of SECTIONS (see Section 2.3) of data 
Total number of patches on the entire dataset 
6 integers indicating the number of each type of 
SECTION in the following order: #NLIF, fLIFT, 

9-word ( 4 ^bytes/word /alphanumeric tit, a 

NSECT sets: 



Record #3: 

I SECT 

NPAT 

NU 

HEAD{15) 

ShSEi^o^PcIStSS^ this section (-HU x HY) 
ti^irr of patches in the "N-line" direction 
.f Hptcfrf* in the "M-l Ine" direction 
l£wonl (dTyfcasJuordi alphanumeric section title 

Record #4: 

THAT (12) 

12 double-precision word trans " 

formtlon Mtrlx {not prt scfttly 

Record 15: P(48) 

• 

• 

Record #(4+NPAT) 

48 double-precision word (8 bytes/wwd) 
PC-patch coefficients In GEwiiRIC form, 
repeated for all NPAT patches on this 
SECTION 


Note that the PRE-PROCESSOR program may be used to generate the PCU input data- 
set for both the on-body data (used in Mode 1, described below) and the flux- 
section data (optionally used In Mode 2, also described below). The PCU data- 
set created by this PRE-PROCESSOR Is relatively inexpensive to create and is 
therefore normally discarded by the potential flow solver. This is true unless 
the potential flow solution is being saved for a Mode 2 case; in the latter 
situation, the PCU dataset is copied into the saved fundamental solution data- 
set created by Mode 1 to ensure that Mode 2 is operating on a consistent geom- 
etry base. 


5.2.2 DF12, Mode 1; FUNDAMENTAL POTENTIAL FLOW SOLVER 

This pro r » Fonts the heert of the syste. In that the fund«ent,l flows are 
generally the most coaplex and expensive part of the solution and form t e 
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basis from which any c<Mbined solutions are obtained. The Input to this pro- 
gran consists of a PCU-dataset (described above. In Section 5.2.1) on unit 1, 
and a formatted card-image dataset which contains some simple control flags on 
unit 5. The format of the control flags dataset is: 


Record #1: 


cc 1-72 TITLE(I), 1*1,18 Alphanumeric run description 

cc 73 MODE A single integer indicating the mode (1) 

Remaining records (as needed) in NAMELIST/Z/ format : 

AREF Reference area (use semi-area if NSYM=1) (default: 1.0) 

B0V2 Reference semispan (default: 1.0) 

CREF Reference chord length (default: 1.0) 

0RIGIN(3) Moment reference center (X,Y, and Z; default: 0.,0.,0.) 

IAUT0W 0 * user-input wake 

1 * automatic trailin g e dg e bisectors 

2 * automatic parallel to x-axis 
(default: 1) 

I COMBO 0 * do not save data for a possible Node 2 omftlnatlon case 

1 * yes, save data 
(default: 0) 

I DEBUG 0 = print standard set of input flags 

1 = print expanded set of input flags (for debugging purposes) 
(default: 0) 

IFUNDP 0 * no fundamental solution printout 

1 * minimum fundamental solution printout 

2 * full fundamental solution printout 
(default: 2) 

IPCY 0 = constant chordwise vorticity distribution 

1 * parabolic chordwise vorticity distribution 
(default: 0) 

IPR132 0 = small print size (164 columns, 89 lines per page) 

1 * large print size (132 columns, 60 lines per page) 

(default: 0) 

IPY 0 = do not save P/V (pressure/velocity) dataset 

1 * save P/V dataset for possible use by ISOPLOT, VECPL0T 
and/or TRACE-ON 
(default: 0) 

IQWIK 0 * do not save a QUIKPL0T”-type output dataset 
1 * save "QWIKPL0T” dataset (see below) 

(default: 0) 


3495H 


45 


QWIKPLOT output dataset is similar to a P/Y dataset in that the data is com- 
pacted into unformatted (binary) form. ISOPLOT, YECPLOT, and TRACE-ON all 
require a P/V dataset to execute; the QWIKPLOT dataset is organized around the 
concept of "strings" of data, where every record was created in the form: 
WRITE(IUNIT)VNAME,N, (Q ( I ) , 1=1 ,N), where VNAME is a double precision alphanu- 
meric string identifier (8 bytes) and Q ( I ) is the string of data. An inhoc'.e 
plotting program (called, not surprisingly, "QWIKPLOT") was written to read 

QWIKPLOT datasets allowing rapid and easy comparisons of results of many CFD 
codes, and/or test data. 

LO 0 = higher-order solution 

1 = lower -order simulation 
(default: 0) 

NSYM 0 = no symmetry 

1 = symmetric about Y=0 plane 
(default: 1) 

ALPHA Freestream angl e-of-attack , degrees (no default; may have up 

to 20 values) 

BETA Freestream angl e-of -yaw (default: 0; must have a value for 

each ALPHA value specified) 

(naming: Catmot use nonzero BETA if MSYM-1.) 

I EXTRA Strip outers of "extra-strips" (see Section 2.6.4), if any; 

these count consecutive strips of LIFT sections only, which, 
as mentioned earlier. If present, must be the first sections 
of the input geometry dataset. 

The 72-panel, si^>le round inlet half-body, (brawn with its image in Fig. 20, 
is supplied along with the program source code as a check case. The coordi- 
nates for this case are shown in Fig. 22. A sample execution of the inter- 
active submit program for a Mode 1 execution of this check case is shown in 
Fig. 23. If the submit program is not used, an input dataset of the form 
shown in Fig. 24(a) must be created by the user. The JCL produced by the 
submit program to execute this Mode 1 check case is shown in Fig. 24(b). 

The output from the MODE 1 execution of this test case is shown in Fig. 25. 

It is basically self-explanatory with the following exceptions: 

P,Q,R the curvature quantities as used in the paraboloidal panel 
definition: z, * PE 2 + 2QEn + Rn 2 

SIGMA the source density value at the control point of the panel 
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VN 

VT 

CP 

CL 

CD 

CSF 

CPITCH 

CROLL 

CYAW 

where 

Vef 

•Vef 

Cref 

q 

L and D 



ETA 


ASTRIP 

SECTCL, 

SECTCD 

CIRCULTH 


the net normal velocity on a panel 


the total velocity Magnitude: V T * + vj + V§ 

the pressure coefficient: 

INCOMPRESSIBLE: Cp = 1 - (Vj/V re f)2 

COMPRESSIBLE: Cp - (p - p re f)/q ref 

the “lift coefficient": L/qA ref 


the "drag coefficient": D/qA^f 

the "sldeforce coeff i cient ■ : Fv/qA^pr 
(Note: F Y is the force in the Y-dirlction) 


the "pitching moment coefficient": My/qA re fC re f 
the "rolling monent coefficient": M^qArefbpef 
the "yawing monent coefficient": M^Aj-eftyef 


* re * Sh00,d * 1»W. 

Megser-lnpul reference span (uhlch should he the seatspu If 


the user-input reference chord length 
the dynamic pressure. pVfef/2 
measured in the lift and drag direction 


“J integrated over the Input panels 

cr ** te<i ,f ss ^ ,) 


on^^e "strip* Y i$ take " aS that of the first control point 


projected (into the X-Y plane) planfom area of a lifting strip 


local strip values of L/qA^ip and D/qAstrip 

of the str "> — 
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5.2.3 DF12, Mode 2: COMBINATION OF FUNDAMENTAL FLOWS 

This program permits the user to combine the fundamental flows (from a MODE 1 
execution) to obtain desired mass flow values (typically within an inlet). The 
combination constants required to obtain the user-defined mass flow values are 
obtained automatically by the program when the user supplies a "FLUX" section 
at the place where the mass flow rates are specified. The cost of generating 
the automatic conbination constants varies linearly with the number of panels 
the user has on his FLUX-section dataset, and therefore may equal or even 
exceed the cost of the Mode 1 solution, although this is typically not the 
case. Optionally, the user may simply input these combination constants 
himself, and thus define his own combination case (perhaps using conbination 
constants obtained from an earlier Mode 2 run). 

Since the fundamental flow solver was designed to handle geometries which con- 
tain lifting leading-edge devices, such as those shown in Figs. 9 and 17, the 
program logic which satisfies the Kutta condition made it necessary to have the 
freestream fundamental flows Include some suction effects as part of the stand- 
ard set of freestream onset fundamental flows. As a result, in order for the 
user to obtain pure freestream onset flows (i.e. without any suction effects), 
a CC = -1 .0 may be used. Note also that up to 5 suction fundamental flows may 
be generated in a Mode 1 case, requiring, therefore, an equal nurt>er of flux- 
setting and/or CC-values to be specified in Mode 2. Furthermore, the nurt>er 
of flux-setting conditions specified may not exceed the nurt>er of FLUX sections 
that are input, although the number of FLUX sections may indeed exceed the 
nurtber of suction solutions available from Mode 1; this latter case is the 
typical one wherein a number of additional FLUX sections are included in order 
to use VECPLOT and/or ISOPLOT to survey the off-body flowfield. 


Since the compressibility correction employed by the present program is the 
Lieblein-Stockman correction which is an "after-the-fact" type of correction 
(unlike, say, the more common Goethert correction which solves a different 
potential -flow problem for each freestream Mach number), multiple Mach number 
results may be obtained from a single Mode 1 set of fundamental solutions. 

Input for Mode 2 consists of the saved fundamental solution dataset fror« a 
Mode 1 case (created when IC0MB0=1), plus a PCU-dataset containing FLUX sec- 
tions (if any), plus a card-image flags dataset, which differs according to 
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whether COHPRS-O or 7. Consider first an Incompressible Node 2 case (i.e. 
C0NPRS=0): 

Record fl : 

cc 1-72: TITLE(I), 1=1, 18 Alphanumeric run description 

cc 73: MODE A single Integer indicating the mode (2) 

(no default) 


Remaining records (as needed) in NAMELIST/Y/ format: 

CONPRS 0 = incompressible flow 

1 = compressible flow (Lieblein-Stockman correction) 
(default = 0) 

IOFF 0 * no off-body points 

1 = off -body points Input on a separate dataset X,Y,Z 3F10. 
(Do not confuse this with FLUX sections which are M x N 
grids of points which produce (N-l) x (N-l) panels; 
off -body points need have no organization into H x N 
grids). 

(default = 0) 

IPR132 0 * small print size (164 col urns, 89 lines per page) 

1 = large print size (132 columns, 60 lines per page) 
(default = 0) 

IFV 0 = do not save on-body P/V dataset 

1 = save on-body P/Y dataset for optional later use by 
ISOPLOT, VECPLOT, and/or TRACE-ON 
(default = 0) 


JPV 0 = do not save FLUX-section P/V dataset 

1 = save FLUX-section P/V dataset for optional later use by 
ISOPLOT, VECPLOT, and/or TRACE-ON 
(default = 0) 

IQWIK 0 = do not save “QWIKPL0T"-type output dataset 

1 = save ■QWIKPL0T"-Lype output dataset 
(default =0) 

(For explanation of format, see Section 5.2.2 on “IQWIK".) 

NCOMB Number of combination cases to be calculated NCOHB values of 

ALPHAC, BETAC, VINF, VREF, etc. must be specified. 

( <20; default: 1) 
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ALPHAC, 
to be 

BET AC 


Requested net -combined- angles of attack and yaw (in degrees) 
dch ieved 

(Note: The program automatically selects appropriate combina- 

tions of the available fundamental flows; however, the user 
cannot request -impossible" combinations, e.g. if all funda- 
mental flows were run with BETA=0, then all BETAC values must 

also be 0. ) , _ ^ . 

(defaults: ALPHAC has none, BETAC defaults to 0) 

Frees tream speed (default: 1.0) 

Reference speed for Cp calculation, 
used for the Mach number correction. 

YINF=0 also, then VREF is set to 1.) 

YC(IC0MB,I) Requested average normal flux velocity, referringto ^ e J.^ 
flux-section, for cognation solution rubber ICOHJ (of NCOro) 
(no default; for CGMPRS-O, either VC or CC must be input for 
each suction fundamental flow generated in Mode 1 ) 

CC(ICCMB.I) Des iota ted amfclnatlon constant far the I-th SUCTION funda- 
mental solution, (default: see YC, above) 


VINF 

VREF 


If C0MPRS=1; then VREF is 
(default: VINF but, if 


A sample execution of the Interactive TSO submit program for a NODE 2 incom- 
pressible case is shown in Fig. 26. If the submit program is not used, an 
input dataset such as that shown in Fig. 27(a) is required to acco^jllsh the 
same program execution. The JCL produced by the submit program to execute 
this Mode 2 check case is shown in Fig. 27(b). 

The output from MODE 2, shown for the 72-panel inlet case in Fig. 28, was 
designed to be self-explanatory and differs significantly from that of MODE 1 

in only two areas: 

1 . The page titled "FLOCMB. FLOW COMBINATION MATRIX DATA" contains the 
details of the automatic computation of the combination constants, which 

are labeled "CC." 

2. For compressible cases (C0MPRS=1), an extra column of the local Mach num- 
ber, labeled MACH, is also shown on the output sheets. 

The input flags for the compressible case (C0MPRS=! ) differ only slightly from 
the incompressible case. In particular, the freestream pressure (total or 
static) and freestream temperature (total or static) must be supplied. In 
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addition, the user Is given the option of specifying either the freestream 
speed (VINF) or the freestream Mach nurt>er (HIHF). The nuatoer of options for 
specifying the flux Is expanded to Include average flux Mach number (MC) or 
average weight flow rate (WC). Finally, the Li eble in-Stockman correction also 
makes use of an incompressible reference velocity (V1BAR) which the user may 
optionally control. 

5.2.4 DF12 Program Limits 

The following program limits and guidelines must be met by the user for ModeV. 


1 . 


Maximum total # 
(this Includes 



2 . 


Maxtmmr total # sections: 100 

(Includes M/MCE sections, etc.) 


3. Maximum total f strips; 300 
(Includes ell sections) 

4 * u^TSeJ onlV LI^an^ T ' S RFY° strips, end -extra* strips (If any)) 

5. Maxlwm # BOLT sections: 5 (see also 11 below) 

6. Maximum # SRFV sections : 5 (see also 11 below) 

7. All LIFT sections (If any) must precede all other sections In the Input 
geometry dataset. 

8. All MAKE sections (if any) must follow all other sections in the input 
geometry dataset. 

9 The order of MAKE sections (if any) must coincide with the order of LIFT 

* sections to which the MAKE sections correspond. 

10. No N-line on any LIFT section may be of zero length. 

11. DBLT and SRFV sections may not both be Input at the same time. 

12. Nonzero BETA cannot be requested If NSYM=1. 


For Mode 2: 

13. Only FLUX sections and/or off-body points may be Input (along with the 
control flags, of course). 

U. Maximum total # FLUX panels plus off -body points : 2000 
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15. The maximum total # FLUX sections : 20 

16. 


17. 


Total # of (VC+MC+WC+CC) conditions specified < (IDBLT + #SRFV) sections 
that were input for Mode 1 . 

"Impossible" flow cob* Inat ions should not be requested, e.a. y jto^ l w 

run with only one angle of attack, then Hode 2 cannot po: 

the Mode 1 fundamental flows to achieve any other angle of attack except 

the one specified in Mode 1. 


was 


5.2.5 Post-Processing Program: YECPLOT 

Input to the velocity vector plotting program, YECPLOT, consists of: (1 ) an 
unformatted P/V (pressure/velocity) dataset (either on-body or off -body, i.e. 
FLUX), and (2) a unit 5 card-image dataset. The "format" of the P/V dataset 
(which Is created automatically for on-body results of MOTE 1 (If IPV*1), and 
either on-body (If IPV-1 ) and/or FLUX sections (If JPV-1) for MOTE 2) Is shown 
in Fig. 30. 


The unit 5 card Image dataset for YECPLOT Is In IWWLIST /IMPUT/ format: 


I DEBUG 0 * (default) normal execution 
1 » generate debug print 

vref Value used to scale velocities before plotting vectors 
(0.0 -► draw all vectors with unit length) 

RYLENG Length of a unit vector in rasters (note that page width, for 
example, is always 4000 rasters) 


NYIEWS Humber of "user-defined" views 
(default: NVIEHS*0) 

KVIEWS(I) Where KVIEWS define up to 10 "standard-views" 

1 * side view 

2 * top 

3 « bottom 

4 3 Inside 

5 * front 

6 * rear 

7 * lower outside front 45° 

8 3 upper outside front 45° 

9 3 lower outside rear 45° 

10 3 upper outside rear 45° 

If KVIEKSd )=0, then all 10 views are drawn 
(default: KVIEWS (1) 3 0) 
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KSECT 


Section ninfcers for which plots will be drawn. Up to 40 
sections can be selected. If no values specified, then all 
sections will be drawn. 


If NVIEWS>0, then NVIEWS additional cards are required: 

PSI(I), THET(I), PHI (I) (3F10.6) 

defining rotation angles (see explanation In TSO submit procedure. Fig. 26) 
for each "user-defined* view. 


A sample execution of the Interactive YECPLOT submit CLIST is shown as part of 
the DF12 Node 2 TSO submit In Fig. 26. The JCl stream that was produced Is 
shown in Fig. 27(c). A sample output of VECPLOT Is shown In Fig. 29. 

5.2.6 Pest-Processing Program: ISOPLOT 


Input to the Isoqram plotting program ISOPLOT, consists of: (1) an unformatted 
P/Y (pressure/velocity ) dataset from NODE 1 or NODE 2, and (2) a unit 5 
card-image dataset containing control flags written In WV&LIST/INPUT/ format: 


I DEBUG 
ISCAL 


IPLOTS 


NVIEWS 


0 * (default) normal execution 

1 = generate debug print 


Scale definition used to set Cp mini 


ISCAL 
1 
2 

3 

4 

5 

(default: 


5) 


-3.0 
-7.0 
-15.0 
No limit 


and Increment values: 
ACp 
0.02 
0.05 
0.10 
0.20 

Automatic 


Plot selection flag: 

0 * generate all plots 

1 * Cp plots only 

2 = delta-star plots only 

3 = skin-friction plots only 


Number of 

defined) 

(default: 


"user-defined" views (up to 9 views can be 
NVIEWS*0) 
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KVIEWS(I) 


KSECT 


Define up to 10 "standard-views" 

1 » side view 

2 * top 

3 = bottom 

4 = inside 

5 = front 

6 = rear 

7 = lower outside front 45° 

8 = upper outside front 45° 

9 = lower outside rear 45° 

10 * upper outside rear 45° 

If KVIEWS(l) * 0, then all 10 views are drawn 
(default: KVIEHS(1)=0) 


Section numbers for which plots will be drawn. Up J° 
40 sections can be selected. If no values specified, 
+h*n ail tprtinns will be drawn* 


If NVIEWS>0, then NVIEMS additional cards are required: 

PSI(I), 1KET(I), mil) (3F10-S1 


defining rotation angles for each "user-defined" view. 

A sample execution of the interactive ISOPLOT submit CLIST is shown is part of 
the DF12 Node 2 TS0 submit in Fig. 26. The JCL stream that was produced is 
shown in Fig. 27(d). A sample output of the ISOPLOT program is shown in Fig. 

29. 


5.2.7 Post-Processing Program: TRACE-ON 

Unlike YECPL0T and IS0PL0T, this is an interactive (on TS0) surface streamline 
calculating program which requires as input a P/V (pressure/ velocity) dataset 
from MODE 1 or MODE 2 and user responses to the interactive questions. In 
using this program, one area which the user must understand is the method of 
telling the program where to "start" streamlines. For the purposes of 
TRACE-ON, the body surface is assumed to consist of a nwd>er of SECTIONS, each 
of which consists of an NU by NV grid of data, where NU and NV represent the 
number of points along N-lines of data, and M-llnes of data, respectively, of 
a given SECTION. Note that both MODE 1 and MODE 2 extrapolate the computed 
potential flow velocity data to the edges of the input SECTIONS. This means 
that there are data values in the P/V dataset than panels printed on the 
MODE 1 and MODE 2 output sheets. For example, say the user inputs a section 
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with N "chordwise" points on each of N H-llnes (see Fig. 31); the ntmber of 
panels produced In the potential flow program 1$ (M-l) x (N-1). But the 
nuaber of data points on the P/Y dataset Is (M+1) x (H+1 ) since the data is 
extrapolated "chordwise" (l.e. In the N-llne direction) at the beginning and 
end of the strip of panels, and "spanwise" to the "Inboard - and "outboard" 
edges of the section. All the data points may therefore be described by 
parametric variables in the N-line and M-line direc- tions; these are referred 
to, herein, by U and Y, respectively. Thus the fist control point of the 
potential flow program is called U=2.» V=2., (not 1 . ,1 . since 1.,1. would 
refer to the corner of the section). The second control point (along the same 
strip of panels In the section) would be U=3., V=2., and so on, as shown in 
the figure. 

A sa^rie execution of the TRACE -Oh program Is shown In Fig. 32, and a plot of 
the calculated streamlines are shown In Fig. 33. 
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Figure 2. Chordwise pressure distribution for the EET configuration 
n = 4.31 c , (a) low order methods, (b) high-order methods 









Figure 3. A plane trapezoidal panel. 



Figure 4. Adjacent panels used in nunerical source density derivatives. 
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Figure 8. 


A general creed surface panel and Its projection In the tan*** 
plane. 
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Figure 9. (c) -Faceted” Solid Rendering of Single-Sl.t Case, 



65 








Figure 9. (d) "Smooth-Shaded" Solid Rendering of Side-View Closeup of 

Single-Slat Case. 
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Figure 1C. 


Comparison Between Ax i symmetric and Three-Dimensional Methods for 
Single Slat Inlet, (a) Pressure Distribution on Cowl and Centerhody 
(b) Pressure Distribution on Slat. 







original page is 
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ero Angles of Attack and Yaw. No Added Suction. 











Figure 13. Flowfleld for Inlet with Single Slat. Combined Solution, 40° Angle of Attack, Zero Yaw 
Fan Face Velocity Twice Freestream. 
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Figure 14. Off-Body Isobars for Inlet with Single Slat. Static Solution. 
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Figure 16. Off-Body Isobars for Inlet with Single Slat. Combined Solution, 
40° Angle of Attack, Zero Yaw, Fan Face Velocity Twice Freestream. 
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Fiaure 18 Comparison of Main Slat Pressure Distributions on Single and uouble 
9 Sla^ Configurations. 40° Angle of Attack, Zero Yaw, Fan Face 

Velocity Twice Freestream. (a) 13° from Top Center, (b) 90 From 
Top Center. 
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Figure 19. 


elocity Vectors for Double Slat Configuration, 
ero Yaw, Fan Face Velocity Twice Freestrean. 
iurface. (b) On Inner Cowl Surface. 


40° Angle of Attack, 
(a) On Outer Cowl 
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Figure 20. Wire- Frame Pictures of the 72-Panel (on the ■'Half-Body") 

Round Inlet. Note the "Doublet Surface" Visible Inside At 
the Rear of (b). 
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Figure 22. Coordinates used fbr the 72-panel simple Inlet check case. 
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H31H6027F: 72-PANEL SIMPLE INLET FUGGED AS "SRFV" SECTION. 
&Z ICOMBO=1 , IPR1 32=1 , ALPHA=0,10 &END 


1 


Figure 24(a). Alternative input dataset for DF12, Mode 1, for the 72-panel 
simple inlet check case. 


//* 

//* PC-PATCH FITTING FOR: H31HG027F 

//* 

//PCPATCH EXEC PGI1=H31IC,PARr1= , l,0,0 , ,REGION=IOOOIC 
//STEPLIB DO DSN=TS0T3DF.Bf!F. LOAD,DISP = SHR 
//FT01F001 DO D5N=T5OT3DF.TE«1PGEOM.DO80985 . T 0921 17 , 

// DISP=< OLD, DELETE) 

/✓FT86F0D1 DD SYSOOT=A 

//FT11F001 DD UNIT=SYSDA,SPACE=CTRK, (30,10)), 

/✓ DCB=(RECFM=VBS.BLKSIZE=19069) 

//FT12F0B1 DD UNIT=SYSDA,SPACE=(TRK,C30.1O)) ,DISP=(NEW,PASS), 

// DCB=(RECFH=VBS.BLKSIZE=19069) 

/* 

//» 

//n 

//* 3-D HIGHER-ORDER LIFTING HEUftAMH SOLUTION. 

t/m CMITH INLET CAPABILITY) 

//* MODE 1. CASCID: G827F 

//» 

//MEUHAHH EXEC P6W=HAIN,REGI0N=4M0K 
//STEPIIB DD DSH=TSOT3DF. DF12 . LOAD,DISP=SHR 
//FT01F8O1 DD DSN=*. PCPATCH. FT12FBB1, 

// DISP=< OLD, DELETE) 

/✓FT02F8O1 DD DSN=TS0T3DF . GB27F . FtlttDSOLN , 

//FT0AFOO1 DD UNIT=SYSDA,DCt=(RECFH=VBS.BLKSIZE=19BG9) , 

//FT05FOO1 DD DSH=TSOT3Df!dF12!dOS 0985.T092117,DISP=COLD, DELETE) 
//F7O6F801 DD SYSOUT=A,DCB=CRECF«=FBA.LRECL=1A9.BLKSIZE=16500) 
✓/FT08F001 DD UHIT=SYSDA.DCB=<RECFn=¥BS.BHCSIZE=19069). 

/✓ SPACE=(TRK,(1,10)) 

//FTO9F0O1 DD UNIT=SYSDA,DCB=<RECFM=VBS.BIKSIZE=19069) , 

// 


SPACE = (TRK,( 1,10 ) ) 


//FTI0F801 DD UNIT=SY$DA,DCB=(RECFP1=VBS,BIKSIZE=19069), 


// 


SPACE=CTRK,C1,10)) 


/ / jrAvc-uAnn*#**/# 

/✓FT11F801 DD UNIT=SYSDA.DCB=(RECFPI=VBS.BLKSIZE=19069), 


// 


SPACE- (TRK, (1,10)) 


//FT12F801 DD UNIT=SY$DA,DCB=(RECFT1=VB5,BLICSIZE=19069) , 


// ii\mi%a»*n// 

//FT13F081 DD UNIT=SYSDA,DCB=(RECFM=VBS,BLKSIZE=19069) . 


SPACE=(TRK,C1.10)) 


// 


SPACE=(TRK.C1,10)) 


// jrNVC-\ 

//FT16F001 DD UNIT=SYSDA,DCB=CRECFH=VBS,BLKSIZE=I9069) , 


// )r out- V I ML 9 V A # J. W S * 

//FTI5F881 DD UNIT=SYSDA,DCB=(RECFI1=VBS,BLKSIZE=19069) , 


SPACE=CTRK,C1,10)) 


// 


SPACE 2 C TRK ,(1,18)) 


// 3rHV»t“v i » s i » i. v j * 

//FT16F0O1 DD UMIT=$YSDA,DCB=(RECFH=VBS.BLKSIZE=19069). 


// 


SPACE = (TRK,( 1 ,10 ) ) 


//FT17F001 DD UNIT =SYSDA , DCB=C RECFPI=VBS, BLKSIZE=19069) , 

// SPACE-(TRK,Cl,10)> 

//FTI8FOO 1 DD UNI T =SYSDA ,DCB = (RECFM=VBS» BLKSIZE 2 19069) , 
// SPACE=(TE.<,(1,10)> 

//FT20FO01 DD DUf'-IY 


Figure 24(b). Mode 1 JCL stream for the 72-panel simple inlet check case. 
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Figure 25. Semple output from the DF12 Node 1 solution for the 72-panel simple Inlet check case. 
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Figure 25. Continued. 
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Figure 25. Continued. 
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Figure 25. Continued. 
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Figure 25. Continued. 
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Figure 25. Concluded. 
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Figure 26. Continued. 
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Figure 26. Concluded. 


INCOMPRESSIBLE COMBINATION RESULTS FOR THE 72-PANEL INLET 
&Y IPR1 32=1 , IPV=1 , NCOMB-2, ALPHAC=0,5, VC-1 .3,1 .5 SEND' 


Figure 27(a). Alternative input dataset for DF12, Mode 2, for the 72-panel 
simple inlet check case. 


PC-PATCH FITTING FOR: TSOT3DF.H31H.G027E.DATA 


MODE 2 , CASEID- GI27F 


//* 

//# 

//* 

ZRE8 IS e g s 527| ii 5 : Sis" 5 |hr”" 1 * ,0K 

//FT06FOOI DD SYSOUT=A 

//FT11F001 DD UNIT=SYSDA.SPACE=<TRIC,C30.10)), 

/✓# 

//* 

//* 

//* 

//» 

//» 

//M 

zsssb SrSJSOT& # Kr?ffi^«~ 

//FT.2F..I dd 8U3[S:&ft3i Si8X; SH ' 

/✓* DISP-5HR 

'/* PV DATASET: 

/ZFT03F001 DD DSN=TS0T3DF.PV .C027F.H0DE2 .OH. 

//K DISP=d»ffl,WTl6),UinT=T50D*,SP*CEi(T1tK,(U,I|),HSE) 

//FTMFMl DO 

"RSSSSi ;; 

Z"“"“ S 

//FT89F001 DD U«T|SYSDA: DCB=<RECFW=VBS,BIRSI2E=1906») . 

//FT10FM1 DD UHTT=SYSDAoDCB=(RECFn=VBS>BLKSIZE = 19069) , 
z/FTllFOOl DD ^NIT^SYSDaIdCB=(RECFH = VBS>BLKSIZE = 19069) o 
ZZFT12F001 DD jniMnN! DCB=(RECFH=VBS.BLRSIZE=19069) . 

^FT13F001 DD S5”g^-JCB?tWCWWS.B«5KE.MM*|. 

Z.FT1AF001 DD jjWTgW g :» B ^CFH=VBS.BLKSIZE=1906„ , 

/ZFT15F001 DD 

z/FTUFOOI DD jjpAgg]f|g£ ’ ^*g*^FM=VBS . BLKSIZE= 1 906 9 ) , 

//FT17F001 DD UHITgYS {{; f CB=(RECFH=V8S .BUCSIZE=19069>. 

//FT1BF001 DD UHIT=SYSDA;dC8=(RECFH=VBS.BIKSIZE=19069) . 

// SPACE=(TRK,C1,I0» 

ZXFT20F001 DD DSN=».PCPATCH.FT12FO0l. 

" DI5P=(0LD, DELETE) 

Figure 27(b). Mode 2 JCL stream for the 72-panel inlet check case. 
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JCL TO VECPLOT A 3-D PRESSURE/VELOCITY FILE 



//* 

//:< 

//* 

//VECPLOT EXEC PC-H = UVECPLT ,REGION = 950K 
//STEPLIB DD DSM = TSOT 3CP . H17 . LOAD. DISP = SHR 

//FT05F001 DD * 

//FT06F001 DD SYS0UT=A,DCB=<RECFM=VA,BLKSIZE=141) 
//FT1SF0 0 1 DD DSN = TSOT3DF.PV.G027F.fiODE2.OM» 

// DI5P=SHR 

//SD9060 DD DSM=ROUTE-DAC.GCf1IF.BON.FL0060. VECPLT . 
// DISP=(NEU.KEEP>. 

// UNIT=T APE16 .LABEL=RETPD=10.DCB=DEN=3 


Figure 27(c). JCL strew to execute the VECPLOT program. 


//* 

//* JCL TO ISOPLOT A 3-D PRESSURE/VELOCITY FILE 

//* 

//ISOPLOT EXEC PGM=UISOPLT.REGIOH=750K 
//STEPLIB DD DSH=TSOT3CP.H17.LOAD,DISP=SHR 

//FT05F001 DD * 

/* 

//FT06F001 DD SYS0UT=A,DCB=(RECFH=VA,BLKSIZE=141> 
//FT18F001 DD D5N=TSOT3DF.PV.G027F.HODE2.ON, 

// 0ISP=SHR 

//SD4060 DD DSN=ROUTE.DAC.GCI1IF.BON.FLOB6B.ISOPLT. 

// DISP=( NEW. KEEP) . 

// UNIT=TAPE16 .LABEL =RETPD=I0.DCB=DEN=3 


Figure 27(d). JCL stream to execute the ISOPLOT program. 
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Figure 28. Sample output from the 72-panel simple Inlet check case. 
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Figure 28. Continued. 
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Figure 28 . Continued. 
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Figure 28. Continued. 
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Figure 28. Continued. 
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Figure 28. Continued. 
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Figure 28. Continued. 


Ift 


I 


Ul I 

44 UJ I 

OC I 
-O I 
«0lU I 
I 
I 

^ I 
30 I 
-> * I 
o I 

I 

I 

• II I 
VO I 
<< I 
O K I 
m ui t 
OC M I 
3 l 
Z * I 


m i 
i 
i 

II I 

0 I 
«< I 

z i 
a. i 

-j i 

«< i 

in i 

z i 

O i 

•>4 t 

V I 

3 I 

-J I 
O I 
m t 

1 

z i 
o I 

1*4 I 

V I 

< I 

*! 

s: 

o I 
u i 
I 

3 i v 
O I Hi 
_f f 

Ik I x 

I M 


OC I IV 
Ul I 

z i at 


Ik I 04*1 
4 I 0 U 
OC I Cu 

o I 03 

oc I om 

04 i m 

< 1 UiUJ 

I «iu 
m I MIL 
4 I M 
-» I 4*0 
O I 0*UJ 
3 • Uk 
O i at < 

a, t^ 


X | 
V I 
-i I 

3 i 
U I 
OC I 

04 | 

O I 


3 * 

4 I 
V I 

u « 


^ 1 

1 

1 

4 1 -4 

• 

IV 

04 i m 

1 

HI 

V 1 * 

1 

OC 1 

Z 1 4 

» 

04 1 

ttf 1 0. 

1 

•4 1 

V 1 1 

1 

NO 1 

O • CM 

1 

»0t 1 

0» i r» 

1 

Oo 1 

• 

1 

m i 

OC 1 UI 

I 

• i 

III 1 x 

1 

IV 

O t V 

1 

HI 

OC 1 

f 

OC 1 

o • oc 

1 

«X I 

1 o 

1 

O I 


0 04 | 

I 


O 1 o- 

1 

0*0 1 

04 | 

1 

UI 1 

z 1 3 

1 

3 

1 01 

1 

.4 

V 1 UI 

I 

< 1 

z 1 oc 

1 

> 1 

4 I 

1 

IV 1 

& 1 z 

1 

V 0 * 1 

C 1 o 

1 

30 i 

O 1 44 

1 

0 . 1 

O 1 V 

1 

Z I 

1 4 

1 

04 | 


QU l 

Ui I 


ORIGINAL 
OF Poe? 




►4 t 

OC I 


l/> I 
■4 I 
I 


O I 
O I 


0» I 
I 


I 

—I I 
U I 
0- I 
O I 
UI I 
4* I 
I 


AIHOOHM 


r*. *r to 

MMNONN 

41^00410 


omomoo 

9>44»NH 

« 0 N«A 444 


nmamnn 

rv Noor^r* 


n«p44«» 


~ ! GSSSSSiS 

«4U4Pt# i « 
Hm 444M « 0 
*44lAA4 4 • M 


I A 

A#«»4p4 I A 
4#04A4 I A 
4N4444 | • 

^o«a«hA I « 

**•♦*• I • 

44AA#4 i M 


AAANAO t 0 

AI^00C«4 | c 


IH04MMO | «* 
IS.NAAN0 I 0 


4*^^^4^0> | O 
I I I I I 


AA0A44 I N 
0ONNO4 I O 
A04NOO | 4k 
AArt^ryf^ i *4 


ONNANO I (M 
II I I I | | 4 


I O 


OOIAAAN I , _ 

| NHHOH0 I 9» II » 
NN4044 | CM 
I • 

Ill II 

** *-• <*4 ^ *n i « 

OONHM4 I o* 
0AA44A | 0k 

0«N4(k« I -4 

0or^0AM I f*. 


it 

0 -A 
H 0 

U « 

M M 

» CM 


0 01 
H 0k 
0 * 
H « 

m m 

m • 
0 N 
0 M 


« M 
4 0k 

0 O 

0 • 
0 o 
4 I 


H fs. 
0 O 
II 0k 
II -4 
N O 
II • 
M CM 
II -4 

I 


H 0* 


II CM 
If • 
II O 
II I 


M 0k 
II 0* 
0 0k 
0 #4 

N I 



OUI 

1 l l 

...... 

• 

•1 * 


oo- 

1 V 1 

NN00NA 

o 

0 O 


zz 

1 1 

1 1 1 


If 


04 04 







1 » 1 

44 CM *0 4- 0k «0 




*® *» 

1 Z 1 


< 

4 


HI HI 

1 1 


*- 

V 


44 



o 

O 

o 

W 

1 O 1 

—4 

V 

k- 

44 

0404 

1 IU 1 





w 

1 0k I 


• 

O 

UI 

UI HI 

1 UI 1 

» 

V 

o 

04 

Ik 

o 

0*0 

i a. i 

U. 

HI 

z 

< 

<< 

1 V 1 

OC 

0) 

o 

a. 

0 4. 

i ►- » 

0k 


o 


P *GE fs 

DUALITY 


117 


Figure 28. Continued. 
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Figure 28. Continued. 
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Figure ZB. Concluded. 
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Figure 30. P/V (Pressure/Velocity) dataset format. 
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2,1 



(l.H+l) 


Figure 31. Extended set of velocity points stored in the P/V dataset: solid 
symbols are standard panel control points, hollow symbols are 
interpolated/extrapolated points. For MxN input nanel defining 
points, (M-l ) x (N-l ) panels are produced, and (M+l ) x (N+l) 
velocity points are calculated and saved. 
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Figure 32. Sample interactive TRACE -ON execution for a streamline on the 72-panel 
simple inlet check case. 
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Figure 33. Two streamlines calculated by TRACE -OH for the 72-panel simple inlet 
checic case. 
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Figure 34. (a) Examples of correct and incorrect input; (b) Plan view of the 

input points on a body divided into sections; (c) Another possible 
division into sections. 






appendix a 

CONSISTENT EXPANSIONS FOR THE POTENTIAL AND VELOCITY INDUCED BY 
A CURVED PANEL AT A POINT IN SPACE 

A "true" panel on a body surface is the ^ The 

boundary curves of t P surface at some central 

def 7!fTeUue paneU^A* paneTcoordinate system Is constructed whose origin 
"•’"I of the trUe ’ 2 or < axis is normal to the tangent plane. 

,s the tan, ency point -d^hos ^ ^ corner points of the 

true* panel Ire projected into the tangent plan. ^ " >m , 

* io T:tz::TZ ir r e .F l « ^ - *. 

II„r ’ts construction ~ ^ 

panel. They are the curves joining the true corner poin 
line projections in the tangent plane. 

IZld in a Maclaurin series, successive terms become small, i.e., 


f (t.n» ■ j„ HT « W ♦ " s Jo Q " (r ~" } 


(A.l) 


tnen 


m < n 


(A.2) 


tipi « IV 1f 

P.SO, the vertical distance c of a point on the true panel is of the order 
0 f the square of the horizontal distance, i.e.. 


= 0(L^ + n^' « / ^ ^ 


(A.S) 


V * Hint (X v z) in space due to a source distribution on the 
The potential at a point U.y.zj in sp 

true panel i r - 


2274H 


A-l 


(A.4) 


♦ ■ fj 

s r 

where 


r ? = (x - O 2 + (y - n) 2 + (z - ;) 2 (A. 51 

and where S is the surface area of the true panel. It is desired to express 
' in terns of a series of integrals over the projected flat panel. Thus it 
■iay be stated that the potential of the true panel is "expanded about" that of 
the flat panel. To illustrate the process more completely, a three-tern expan- 
sion is derived, although only the first two terms are ultimately retained in 
the present method. 

In panel coordinates the equation c = f(£,n) of the surface of the true 
panel may be expanded in a Haclaurin series in the form 

t = [P£ 2 + 2QSn + Rn 2 ] + [T 3() S 3 + T 2T S 2 n + T^n 2 + T^r, 3 ] + ... (A.6) 

- ^ ^ ^3 + • • * 

There are no constant or linear terms in (A.6) because the origin is at the 
tangency point. All coefficients in (A.6) are constants proportional to 
derivatives of c at the origin. The coefficients P, Q and R, which are the 
only ones actually used in the present method, are the second derivatives. 

They are referred to below as the surface curvatures to which they are closely 
related. The quantity c 2 is second order in S and n and thus in panel 
dimension t, and s is third order. 

The equation of the true panel nay be written 

F(S,n,0 = c - ' s 3 U,n) - ... = o (A. 7) 

Then, taking the gradient gives 

grad F = £ - ^ + ...)i - + A 3n + •••)! ( A.8) 

v/here subscripts F and n denote partial derivatives and f, j, fT are the 
unit vectors of the panel coordinate system. The vector grad F is normal to 
the panel at any point. The unit normal at any point is 


A- 2 


2274H 


(A. 9) 


t _ grad F 
n " | grad FT 


so the ^-component of the unit normal is 


"c = I grad FT 


(A. 10) 


/ 1 t (C 2c t ; 3i; * -) Z ♦ C 3n + —l 2 


This can be expanded in the form 


\ - ' + ^3? + -> 2 + U 2n* C 3n * ••• )2j 

* I lk 25 + c 35 + '•• )2 + U 2n * *-3n * ••• |2j2 + 


(A. 11 ) 


From (A. 6 ) it is clear that c 2 ^ and ; 2n are first order, c 3 ^ and c 3ti are 
second order, etc., so the leading terms from the first square bracket of 
(A. 11) are second order, and the leading terms in the second square bracket 
are fourth order. Thus, 


i 1 , ,2 2. 

n C " 1 ' 2 + W 


(A. 12) 


is a valid three-term expansion with second (linear) term zero. More precisely. 


nj. = 1 - \ L(?PS + 2Qn) 2 + (2Q£ + 2Rn) 2 J 

= 1 - 21 (P 2 + Q 2 )^ + 2(PQ + QRKn + (Q 2 + R 2 )n 2 J 


(A. 13) 


n^ = 1 - 24- 2 

where is second order. The elementary surface area dS on the true panel 
is related to the elementary area dA = d£dn in the tangent plane by 

n^dS = dA (A. 14) 

d$ = (1 + 2* 2 )dWn (A. 15) 

The source density is strictly a function of surface distances along the panel. 
However, it can be shown that there is no difference between surface distances 
and £ and n through second order. Thus it suffices to define o in terms 
of r, and n in the form 


A- 3 


2274H 


o 


= a 0 + (o % K + 0 y n) + (o xx ^ 2 + 2o xy ^n + y 2 ) + ... 

= a + a, + o- (A. To) 

o I Z 

where o n is n-th order in £ and n and thus in panel dimension t. The coef- 
ficients in (A. 16) are constants proportional to derivatives of o at the 
origin. 


Thus, to three terms 


odS = (o Q + o-j + 02 ) (1 + 

= (o + o, + o' )d^dn 
0 1 L 


(A. 17) 


where 

o£ = 02 + 24»2° 0 (A. 18) 

All of the above expansions are independent of the location of the point 
(x,y,z). 


It remains to expand 1/r, which of course, does depend on x,y,z. It is neces 
sary to differentiate three ranges of value of r: 


(a) r » t 

(b) r = 0(t) 


r « t 
(c) 

r = Olt) 


for all 5,n 
for all 5,n 

for some £,n 
for other 


(A. 19) 


In other words, the ranges amount to the situations where the distance of the 
point (x,y,z) from the origin of panel coordinates is, respectively large, of 
the order of, and small compared to the dimensions of the panel. It turns out 
tnat the range (b) where r is of the order of the panel dimensions is the 
essential one. In the far-field, range (a), 1/r is expanded in negative powers 
of r Q * (x 2 + y 2 + z 2 )^ 2 . The resulting expansion differs in no important way 


2274H 


A-4 


(b) shows that while certain valid for range 

— r morels, the expansion - ^ could be elimi- 

(a), although It is conserve ,n that some r t0 assume 

nated. For points Cose to the panel rang Vi , ope «ith respect 

** ’T, ^"^tr^rr^-magnuude ana, y sis oF 

to the tangent plane Under t ^ ^ of , panel method 

nrr r;r~ ; “r;i:r,r:7r»r;,.,, 
ErBsrsrrr- 


The distance r can be written 


r 2 = ( x - l) Z + (y - + z 2 - 2zc + S 

= - 2z£ + C 


(A. 20) 


r f is the distance between (xjr.rt and the point «.«.« » the flat 
panel (Fig. 3). Thus 

. ’n 1 -2*t + ^ 4 ^5- + -••! 

i . i I * tt - n - 7 — 3 + s 

r ” r f , TT, Z r f r f r f 

f / 1 + (-2zc + C )/ r f 




(A. 2 1 ) 


Note that (z/r ) is of order unity, and, since r is 0(t), the quant y ^ f 
s also 0(1 k L. 0(0. K -e is .ade of (A.6), the above bec.es 


A-5 


2274H 


(A. 22 ) 


!.L.n + f_ 

r r/ 1 r f 


Co + C, + - 


3z 2 h.*2*h* 


* ' 2 ' '’3— - + (|Vl )( 

r. 


2 

) ] 


f 'f 

from which the desired three- term expansion is obtained in the form 




2 

1 * "2 
7>"? J} 


(A. 23 ) 


or, for abbreviation. 


1 . 1 (1 * c. ♦ C 2 ) 


(A. 24) 


where c is of order n in ? or n and thus in t. Now this is multiplied by 
n 

(A. 17) to obtain the expansion 


£ dS = -J- (o Q + + o')(l + c-j + c 2 )dA 


= r [o o + °1 + °2 

r 


(A. 25) 


+ c 1 o 0 + c 1 o 1 + c^ 


+ c„o_ + c„o, + c„o;] 


'2 o 2 U 1 2 2 

The square bracket in (A.25) must be reduced to a three-term expansion using 
the facts that 


and 


o » o. >> o« 
o I ^ 


1 » c, » c 


1 >> c 2 


(A. 26) 
(A. 27) 


The leading (lowest order) term is clearly o Q because all terms are small com- 
pared to it. Possible members of the second term are o 1 and c because all 
remaining terms are small compared to one or the other. Thus, both of these 
are retained in the second term, because neither can be guaranteed to be small 
compared to the other in all cases. The question then arises as to which of 
the remaining six terms of the square bracket of (A. 24) need to be retained 
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the third term of the expansion. Obviously, if any of these six is small com- 
pared to any other, it may be discarded. This eliminates all but o^, c^o-j and 
c 2°o' Thus, the three-term expansion of the integrand of (A. 4) has the form 

F dS s 77 [o o * (c l°o 4 V 4 (c 2°o 4 Vl 4 "j’ 1 


(A. 28) 


when the abbreviations of (A. 16) and (A. 24) are replaced by their actual 
expressions, the three-term expansion of (A. 4) is 


[• 


2 2 

Co o o Co 1 Co Co 

r [j (z -3 + g z "§ • ] + P 2 "T {o x*> + a y n)dA 


(A. 29) 


+ f J f- + 2o xy ^n + o yy n 2 + 2 (P 2 + Q 2 )? 2 + 4(PQ + QR)£n 

A f 

+ 2(Q 2 + R 2 )n 2 ] dA 

where the integrals are over the projected flat panel. Defining the integrals, 

jn n 

I nn - fj ^-dA (A. 30) 

mnp J J r p 


the three-term expansion of the potential can be written 


4> = 4» (0) o 0 + [$ (c l 0 


.dx). 


a„ + '"'o v + 4>' ,J 'oJ 


.Oy). 


+ [«f (20) o 0 + $ (2x) a x + <t> (2y) o y + (J) 


y J 

(2xx) , 
xx 


„(2xy) . 


♦ 2<p' t '' J 'a* + ] 


,(2yy) , 


(A. 31) 


xy 


yy 


where 


o' = a + 2(P 2 + Q 2 ) 
xx xx v 

c' y = o xy + 2 (PQ + QR) (A. 32) 

°yy ' °yy 4 2 "> 2 4 r2 » 
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(A. 33) 


The individual potentials in (A. 30) are 


*<°> . i 


001 


4>^ = z[PI 2 03 + 2( 5 I 113 + RI 023^ 


* (1<) - I 
♦ " 1 101 


JlY> - I 

* ' ^il 


(A. 34) 


(A. 35) 


(|)^ 20 ^ = z[T 30 I 303 + T 21 I 123 + T i2 T 123 + T 03 I 033^ 

* I z 2 [P 2 I 405 * «QI 315 * (2PR * V)I 225 ♦ WI,35 * r2 W 

- 7 [p2l 403 + 4PQI 313 + (2PR * 4q2,I 223 * 4<1RI 133 + r2i 043 ] 

(A. 36) 

<^ 2 ** = z[PI 30 3 + + RI 123^ 

(A. 37) 

*^ 2X * = z f pi 213 + 2C * I 123 + RI 033^ 


(^(2xx) = j 


201 


J2xy) 


111 


,(2yy) 


021 


(A. 38) 


The first term of (A.31), <|>^o 0 , corresponds to a flat panel with a constant 
source density. This is, of course, the term used in the first-order method. 
The second term of (A.31) contains the second derivatives P, Q and R, of the 
surface shape but no higher derivatives and first derivatives of the source 
density but no higher derivatives. Thus the second term of (A.31) corresponds 
to a paraboloidal panel shape with a linearly varying source density. The 
third term of (A.31) contains all the preceding quantities and also the third 
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derivatives of the surface shape and the second derivatives of the source den- 
sity; a cubic panel with a quadratic source density. 

The equations above Illustrate the fact that succeeding tenns in the 'T*™ 
for the potential of a panel increase rapidly in complexity A single 

Term is followed by a second term Individ- 

its own integral of the form, Eq. (A~>0). The tn.ra 

u.l parts which together involve 17 different ^ ° the ^ 1rd 

(A. 30). The great increase in complexity associated with retaining t 

term of (A.31) appears to be unjustified at this ti«- ^ the 

higher-order -hod accounts for *. source ^ 

zrs:z r— - — * — iRets - 17 

and 18). 
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APPENDIX B 

GENERATION OF PANEL GEOMETRIC QUANTITIES BY MEANS OF BICUBIC SPLINES 


Very elaborate geometry fitting procedures based on parametric bicubic splines 
have been developed at Douglas Aircraft Company over many years. A description 
of this technique is beyond the scope of the present report. A survey is con- 
tained in Ref. 19. In the present application the method is considered a 
"black box," although several minor changes had to be made. 

The points defining the body are input in the usual way. Each panel is fitted 
by a bicubic surface in terms of two parameters, u and v, that vary from 0 to l 
over the panel. (The panel is the unit square in parameter space.) This per- 
mits the well-known procedures of Ref. 20 to be used as follows. 


Let a point, (x,y,z), of the panel be represented as a vector 

-*• -*>. -*■ . 
x * xi + yj + zk 

The parametric cubic fit then yields 

x = x(u,v) 

These expressions may be differentiated analytically to give 

-*■ + -+ ■* 

X, X, X, X, X 

U V uu uv vv 


(B.l) 


(B.2) 


(B.3) 


as functions of u and v. The vectors x u and x y + are tangent to the curves, 
v = constant and u = constant, respectively, and thus lie in the surface 
although they are not perpendicular. 


The point corresponding to u = v = 1/2 is in the "center" of the panel in some 
sense. It is selected as the control point and origin of coordinates of the 
flat projected panel. The derivatives of Eq. (B.3) are evaluated there, and in 
all that follows x and its derivatives are assumed to be those at u = v * 1/2. 


The unit normal vector to the panel, which is also the unit vector along the 
5 axis of panel coordinates is 
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(B.4) 


n = k ^ - + 
e — 


*u x x v 


X U x x v 


where the sign is selected tc give an outward normal. The unit vector along 
the z, axis of panel coordinates is taken tangent to the v = constant curve 
which nearly parallels the N-lines, 


T = “ 

• m 

Thus the unit vector along the n axis of panel coordinates is 

L = It, x t_ 


(B.5) 


(B.6) 


The components of the three unit vectors thus obtained comprise the transfor- 
mation matrix. 

Now define 


h * u - 1/2, 


k * v - 1/2 


(B.7) 


and consider the Maclaurin series for E, n and c in terms of h and k. 
They have the form 

z, - Ah + Bk + (second order) 


n = Ch + Dk + (second order) 


(B.8) 


x, - l/2(eh 2 + 2fhk + gk 2 ) + (third order) 


There are not constant terms in (B.8), because the origin of panel coordinates 
corresponds to h = k " 0. Furthermore, since the plane is tangent to the 

surface at the origin (1c e is the normal vector), the series for x, has no 
linear terms. Reference 20 gives the coefficients of Eqs. (B.8) as 


A = x *i 
* x u e 


B = x. 


(B.9) 


2274H 


D * Xy * ^ 


c 



e - * uu • " f = *uv * " g = * vv 
The first two of Eos. (B.8) may be inverted to give 

h = aE + bn + (second order) 
k = cE + dn + (second order) 


n 


(B.10) 


(B.ll) 


where 


a 




(B.12) 


Equation (B.12) may be inserted into the third equation of (B.8) to give the 
desired form 




PC 2 + 2QEn + Rtf 


(B.13) 


The result is 

P = 1/2 [ea 2 + 2fac + gc 2 ] 

Q = l/2[eab + f(ad + be) + ged] (B. 14) 

R = l/2[eb + 2fbd + gd 2 ] 

For generality c has been included in Eq. (B.14), but in the present applica- 
tion it is zero, which simplifies (B.14). 

It remains to compute corner points in panel coordinates. The four input 
points bounding the panel are transformed into panel coordinates to obtain 
(£*, t£, j£)» k = 1, 2, 3, 4. They are projected into the plane by simply 
ignoring Next the side between points 1 and 2 is rotated to make n-j s hg* 
The midpoint and length of the side are, respectively, 
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1 ★ * 

K = o (r., + r, 9 ). 


• I W » 

n = 2 + 


/, * *.? , * * ; 
d = / ( C 1 - K z > + (n 1 - n ? ) 

Then the final corner point coordinates are 


(B. 15 ) 


n 2 = n 


«1 = 1 - I 


(B.16) 


^2 + ? 

A similar calculation is performed for the side between the points 3 and 4. 


It should be noted that the underlying parametric cubic geometry routine uses 
the surrounding input points to generate the fit to a panel. The routine con- 
siders only points on the same section, and thus slightly different results can 
be obtained depending on how the body is sectioned. For fitting purposes, the 
wake is considered a separate section, so that the routine does not try to fit 
around the trailing edge. On the semi-infinite last-wake panel the derivatives 
P and Q are set equal to zero, so that the panel has straight generators in 
the stream direction, but R, the spanwise second derivative, may be nonzero. 
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appendix c 

AREA MOMENTS OF A PANEL 


Th e normalized mm*nts of the area of the tangent pane, are reared. These 
are defined by 


I = - Ltpt \\ E n n ra dEdn 
1 nm ^n-HTH-z 


(C.l) 


where the region of integration is the area of the pane,. For «anp,e^>oo 
. % t. t 4 I are the moments of inertia or second moments, 

is the area t Un* t l 02 arp two first- 

The order of a moment is the sum of its subscripts n * m. There are two 

r e moments, three second-order, four third-order and five o^ -or^ . 

The present method uses up through fourth order. The moments are ca,cu,ated 
by a straightforward but rather lengthy set of formulas. 

First, normalize the corner point coordinates by the mat, mum diagonal, 

\ - V l > \ * V 1 * k = 1, 2, 3, 4 (C.2) 

now the normalized moment may be defined in terms of certain auxiliary 
functions 

•m+1 f jn+l f n+ ^l 
+ .H3 (S4 " ^3 ,J 

The auxiliary function is as follows: 


(C.3) 


If |r»32l > 1 : 

(32) 
l nm 


JZ) _ 1 - r r n+ 1 JJ n+ 1 ] 

a “ ( m + 1 ) ( n + i ) 


n+lMa+l-iZ 
3 


1 _J_ rr n+2 n n l 2 

( n + 1 i ( n + zT 0132 


n 


+ Tn + Din + LS " 3 


1_ L ^n+3-n-l 3 2 
32 


(C.4) 


n(n - 1) 1 rr n+ 4«n»-2-|2 

(n + TTTn + ? ) ( n + ^)(n + TT 3^ U 

J £_ 


3 


n(m - 1 ) (m - 2) 1 r-n+5^-3^ 

+ (n + 1 ) (n + 2)(n + TRn - 4)(n TTJ 4 ^ n J 3 

™32 


n(n - l)(m - 2)(n - 3) 


(n + 1) (n + 2 ) ( n + 3)(n + 4)(n + 5){n + 6) ^ 


1 r »n+6»n-4-,2 

~T~ U n ] 7 


nn (n + T ) (m +~ZJ ^2 U n J 3 


if i 11132 ! 1 1 
t (32) = 

n „2 r ^n-l»rH-3-|2 

' (m + Dim + 2Jlm TT[ "32 ^ 11 J 3 

n(n - 1) „3 r .n-2yH4-,2 

+ 3){m + '*T ”32 U 11 J 3 


+ Ti + 1 Hm + 2Tn» 

n(n - 1 )(n - 2) 


(m + 1) (m + 2) (n + 3)(m 


1 4 r »n-3»nw5-,3 

+ 4 M n . ~TT "32 t!: n 3 


(C.5) 


n(n - l)(n - 2)(n - 3) ^ 

+ (m + l)(m + £)(m ♦ 3){n ♦ 4 )(m + 5)(ri +TT m 32 15 J 3 


where the bracketed symbols are defined by 

tsW 3 = ^ ‘ c - e > 

(The superscripts in the above equations denote simple powers of the quantities 
except for the bracketed double superscript (32), which denotes the side^of the 
quadrilateral.) It is clear from the above that the calculation of 1^ 
requires n + 2 terns of Eq. (C.4) or n + 1 terms of Eq. (C.5). The calculation 
is simply terminated at this number of terns. The auxiliary function I n|n 
is obtained from the above by an obvious substitution of subscripts. 
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APPENDIX D 

NEAR-FIELD SOURCE FORMULAS 


If r /t < p 2 , the near field formulas are used to compute induced velocities. 
The calculation starts with the element coordinates x, y, z of the field point 
and the geometric quantities associated with the element that are discussed in 

Section 2.3. 


Preliminary quantities to be calculated are: 


r [; = /(x - K v ) 2 + (y - \) 2 + * > 
x - 

a k = ~T. * 


6 -Hi 

h ' r. 


z 

Y k = rT 


k = 1,2,3, 4 

(D.l) 

k = 1,2, 3, 4 


pr 3Z) = ™n [zZ + (y ' n k )2] " (x " ** ){y " T * ) * 

Pk 41) * m 41 [zZ + (y ' \ )Z ^ " {x " ^k )(y ' 


k = 3 or 2 

(D.2) 

k = 4 or 1 


”he basic functions are 


. (mn) _ , n n jam m,n consecutive, i.e., mn - 12, 32, 34 or 41 

l lug r + r + d„„ f 


m n mn 


(D.3) 


and 


(32) 

T, (32) = tan -1 [— — — ]» 


zr. 


_(41 ) 

t ( 41) . tan -l ], 

K zr 


k = 3 or 2 


k = 4 or 1 


(D.4) 


(32) 


Also needed are derivatives of the T's and L's. The derivatives of T k are 


3T 


(32) 


3x 


z,r k 6 k * '’k 321 ^ 1 

^171 
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(oo\ x 2 (32) c 1 

3T k *" zt(2m 32 B k - a k )r k - P k \1 


3y 


3T 


(32) 


XJ2T 


2m 32 z Z r k - p[ 32) (^_zV 


3z 


T&T 


k = 3 or 


(D.5 ) 


There 


D« 32 ' ■ z 2 r 2 ♦ [p> 32) f 

- t (4D 

is an analogous set of formulas for the derivatives of T fe • 


The derivatives of L^ nn ^ are 

(mn) a. (mn) n fR + c \ 

“~ 5 x - D nn (a n + a n K ' "" ( m " ’ 

2d mn 

Dn " ~K * r n )l ' d mn 
mn « 12, 32, 34, 41 

The flat-panel constant-source velocities are 

v <0) _ _^_ l (32) + ’ l (4,) 

S 41 


Jmn) , . 

°L = n v + y l. 
— 5 T~ u nn 1 ’m T n'’ 

(D.6) 


32 


v (0) . _ l (12) + l (34) + ^32 l (32) . L (4D 

V y ^32 S 41 


,(01 = . t <32> * t <32) „ ,««) . T < 41 ' 

2 L ^ 


(D.7) 


Referring again to Appendix A, it can be seen that the integrals I mp of Eg. 
(A. 30) are source potentials if p = 1 and, when multiplied by z, are dipo 
potentials if p - 3. Specifically if represents the potential of a 
dipole distribution v - s’V on the panel, then 


<b = Z T 

y mn 


(D.8) 


mn3 
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It turns out that the higher-order source terns for a panel are expressible in 
terns of derivatives of the dipole potentials. Eg. (D.8), and the derivatives 
of the source velocities, Eq. (0.7). 


Only the derivatives of and V y are needed (since V z = <t>Qg, its 
are exactly a potential derivative). The derivatives of V x and V, 


3 V ( 0 ) - 

i 

3 L ( 32 ) 

X 2 - 

a (4l > 

3x 


3x 

^7 

3 x 

1 

S X 

i 

a (32 > 

+ 1 

ai< 4,) 

ay 

^32 

ay 

^41 

ay 

3V (0) 

X 

1 

3L (32 > 

1 

X -- 

S l< 4 '> 

8 z 

S 32 

3 Z 

S 41 

3 Z 


aV<°> 

3L (,2 > 

+ 

3L (34) n 32 

a <32 > 

_ "41 

a L (4,) 

3x 

ax 

ax 5^ 

3x 

^41 

ax 

av<°> 

3L (1Z > 

x 

a* 341 “32 

a< 32 > 

m 41 

aL< 4 " 

ay 

ay 

* 

3y + ^ 

ay 

ay 

«»> 

a" 21 

X 

al _(34) , m 32 

ai< 32 > 

"41 

a L < 4 '> 

3Z 

3z 

T 

3z + 

3z 

az 


Now the potential derivatives are as follows. 


9 *oo 

aTj 321 

X 

ar' 32 ’ 

+ 

aij 41) 


~1T = " 

8x ' 

¥ 

Sx 


Sx 

3x 

a4»oo 

ar< 32 > 

x 

S t« 32 > 

x 

aT< 4 '> 

aT< 41 > 

ay 

ay 

T 

ay 

T 

ay 

ay 

a> 

0 

1 

ar< 32 ' 

x 

ar< 32 ' 

x 

8T< 4 " 

ar' 4 " 

3 Z 

3z 

T 

az 

T 

az 

dz 


2274H 


0-3 


derivatives 

are 


(D.9) 


(D. 10) 


34>, 


01 


3x 


= - z 


3 V i 0) *00 

-— + y — — 


3x 


3x 


34> 0 i 


3V 


( 0 ) 


ay 


= - z 


y . .. ~ v oo 


a«t>r 


ay 


+ y 


3 + V z ( source ) 


34 


01 


3z 


= - z 


av (0) 

av 34> nr 

3 7 -- + y -5^ - V (source) 


34> 


10 


3x 

34> 


= - z 


10 


ay 

34> ]0 
1 z~ 


= - z 


3V[ 0) H 00 

-*r~ * * — + V z (source) 

aV (°) aA 

3V x + **00 

* X 


ay 

( 0 ) 


3V x 

"5z 


ay 

94>r 


- z — + x - V x ( source) 


Now define 


'11 


r 3 - r r i - r 4 m 32 


2 — + — ~rr ~ + ;y (x ' " 32 * ' b 32 ,L 


(32) 


32 


S 41 S 32 


m 


41 


- rr (x - m 4i y ‘ b 4i )L 


(41) 


41 


3x 


ay 


“3 " a 2 

”T 2 

k m 32 

h 7T 

L (32) 

«!>32 

+ 73 " 

(x - 

■ m 32y - 

b 32 ) 

3L (32) 

3x 

j2 

^32 


" > 32 





a l - a 4 

m 41 

. (41) 

ra 41 

(x - 

- <n 41 y - 


3L (41) 

72 " ‘ 

■7T 

L 

' S 3 " 

b 41 ) 

ax 

S 41 

S 41 


S 41 





63 - e 2 

m 32 
‘ ^ " 

L (32 ) 

tn 32 

i 7 " 

(x - 

■ m 32y ■ 

b 32 ) 

3L (32 > 

3 y 

j2 

" > 32 


*32 





h ~ h 

7 ? — ‘ 

m 41 

L (41 ) 

m 41 
‘ 73" 

(x - 

■ m 4iy - 

b 41 ) 

3l< 4 ” 

~w~ 

S 41 

S 41 


S 41 






( 0 . 11 ) 


(D. 12) 


(D.13) 


(D.14) 
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aL (32) 
' ' 3z 


9J I1 y 3 ' y 2 

~zr -77 

32 


m 


32 


32 


(x - ta 3z y - b 32 ) 


Y 1 • Y 4 


41 


m 


7T (x - m 4] y - b 41 ) 3L 
41 


(41) 


3z 


Using the above 

(Q) _ 


- - c -4r + x ^ * y - *» t - «i 0, 3 


f (Q) _ 




11 


9 4>i 


- - 9 ir + * -sr + y 

3J 


9<J> 


10 


3<l bo ,„(0) 


3y~ T A ~37" t y ~zy~ “ xy ~W~ ~ zV x ] 


V (Q) _ r , J 11 , „ a <bi . .. 34, 10 3<t> oo , 

V z - [z Tf + x “3z“ + y ~3z" ■ xy -~5z" + J ll ] 


Also define 


r 2 ” r 3 


r 4 ' r l 1 


H 02 = ra 32 72 + m 41 12 + IT {x ■ m 32 y “ b 32 )L 


(32) 


32 


S 41 S 32 


’ 7T (x ‘ m 4l y ’ b 41 )L(4]) 
b 41 


aH 02 m 32 


~95T = 72~ (a 2 " a 3 l + TT 
b 32 s 32 


7T vu 4 - V - rr 

h! S 41 


!H02 _ ^2 m 32 , ( 

3y " <-2 [ *2 " B 3 “ 3“ L 
b 32 s 32 

n»4i m 41 

+ 7 T ^4 “ e l } + T“ 
41 *41 


a ) + 1 L (32) + U I ro 32 y ~ b 32 } 3L 


(32) 


? 

5 32 

3x 

(x - m 41 y - b 41 ) 

3L (41 > 

$3 

5 41 

9x 

(x - "32 y - b 32> 3L 

(32) 

c3 

b 32 

ay 

(x - n 41 y - b 41 ) 

3L (41 > 

+ s 3 

*41 

ay 
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(D.15) 


(D. 16) 


(D.17) 


\ 


9H 02 _ 32 ( % 

ST'-jr (y 2 y 3 } 

b 32 


m 41 . . 

+ TT (y 4 - Y l } 
b 41 


(x - ro 32 y - b 32 ) 9 j_(32) 

73 


9z 


32 


(x - m 41 y - b 41 ) 3L (41) 

„1 3Z 


41 


Using the above 

V CR) 

X 

„(R) 


= - [z 


aH 


02 


ax 

aH, 


_ 3<1> oi 
+ ^ y ax 


[z^ + 2y 
3H, 


3y 


(y 2 W 2 ) ^55 ] 




34> r 


. (y 2 * z 2, ^00 . 2zV lO)j 


,( 0 ) 

z 


Finally, define 


i - u _ 

J 02 “ H 02 zv z 


1 - Oj(0) ^ U 

J 20 " * H 


02 


3J 20 _ v (0) ^02 
ax " x ax ’ 


^20 . _ v (0) _ ^02 
ay y 3y 


aJ 20 v {0) 3H 02 
az ' z " az * 


where 


IX ~ C2 ’ ^ 2<y -J*V 32 > 


- (y - n 3 )L 


(34) 


(x - V - m 41 (y - n 4 ) l(41) _ ^(Q) 


41 


Using the above 

„(?) 


aJ 2Q W 10 2 8<l> oo 

~Tix“ + ?x -ar- ' x “ax“ 


2zV* 0) l 


(D. 18) 


(D. 19) 


(D.20) 


(D.21 ) 
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% 


/ (p) - - 
y 

„(p) . _ 


[z 

[z 


sJ 


20 


d 4 > 


ay 

aj 


+ 2x 


10 

ay 


20 


34 > 


10 


2 H 00 
x ~W~ 1 

2 8<t> 00 


+ 2x - x- - 5 J- - -20 


"5z“ 

V (lx) = xvi 0) - J 


+ Jon ] 


20 


v (lx) . xV (0) 


- J, 


y 'y 


y(lx) _ x y(°) _ z v 


( 0 ) 

X 


»(iy) 


yv* 0) - J 


v (iy) = yV 


( 0 ) 


- o 


li 

02 

( 0 ) 


u(iy) = - zv 

v z y z y 


(D.22) 


(D.22) 


(D.23) 
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APPENDIX E 

INTERMEDIATE-FIELD SOURCE FORMULAS 


If 

Pi > r Q /t > p 2 


the intermediate-field formulas are used. 
First define direction cosines 



V 

b = 

r o 

y 

It 

cl" 


(E.2) 

Next define certain " 

derivative 

functions" as 

follows: 




First Order: 








C 

X 

II 

• 

Q 

w 


Uy = - 

-B, 


u z = ■ 

-y 

(E.3) 

Second Order: 








u xx ■ *** ■ 

1. 

u xy = 

3aB, 


Uyy * 

3B 2 - 1 

(E.4) 

U XZ ■ ^ 


u yz = 

3By, 


u zz = 

3y 2 - 1 


Third Order: 








u = 3a(3 

XXX 

- 5a 2 ), 

u xxy 

= 3B(1 

- 5a 4 ), 

u xxz 

= 3yd “ 

5a Z ) 

u = 3a( 1 

xyy 

- 5g 2 ) , 

u xyz 

= -15aBy, 

u xzz 

= 3a(l - 

5y 2 ) 

(E .5) 

u = 3g(3 

yyy 

- 5D 2 ), 

u yyz 

= * i 1 
' 1 ' 1 

- 5B 2 ), 

u yzz 

= 36(1 - 

5y 2 ) 


Fourth Order: 

u = 9 - 90a 2 + 105a* 
xxxx 

u = 15aB(7a 2 - 3) 
xxxy 
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u = 15ay(7c/ - 3) 
xxxz 1 

u = 3 - 15(ct 2 + p 2 ) + 105a 2 6 Z 
xxyy 


u xxyz * 156y,7 “ ' 11 

"xxzz * 3 ‘ 15,1,2 * ■ ,2) * l05a2y2 


(E.6) 


“xyyy ' l5aB(7B ' 31 


“xyyz = 1W7B • 11 

“xyzz ■ ,5 “ B,7y2 - 11 
u = 9 - 90 p 2 + 105p 4 

yyyy 


U yyy 2 = 15P Y (7P" ~ 3) 

u yyzz = 3 - 15(p 2 + Y 2 ) + 105p 2 Y 2 


The source velocity components are: 

2 2 
v i 0) = T {-^x + ^ r - ^ ^ 1 1 O u xx * 2 ^“xy-* " 2 ^20 u xxx + 2I 11 u xxy 

m a n 


+ I 02 u xyy^ 

2 2 
Vy°^ = Tf Moo u y + ^F~^ I 10 u xy + ^“yy^ ' 2 *7“* ^20 u xxy + 2I ll u xyy 


•y " z 1 *oo M y ' 'r ,L 'l0* I xy ’ ‘Oryy- - . Q 
0 


+ *02 u yyy^ 


(E.7) 


y(°) = t^ { _, + (£-)[I 10 u xz + IgjU 3 “ \ ^ ^20 u xxz + 21 ll u xyz 

r Q ° o 

+ Vyyz^ 
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Far-Field 


1st Order 


2nd Order 


*+ z 

V X Q = ^11 u xz “ ^r - ^*21 u xxz + *12 u xyz-^ + 7 ^31 u xxxz + 2 *22 u xxyz + *13 u xyyz-^ 


V (Q) ~ t , * 
v y ^ {I ll u yz 
r o 


V Z Q) = "T {I 11 U 2 


t 1 t 2 

(^)[l2] u X y Z + J^yyz-l + 7 ^31 u xxyz 4 2I 22 u xyyz + ^“yyyz-^ 
° ° (E.8) 

t It 2 

( r~ )[I 21 u xzz + I 12 u yzz ] + 7 { r _) [I 31 u xxzz + 2I 22 u xyzz + ^“yyzz^ 
0 0 


4 2 

V x T ^02 u xz " ( ?“ ,[I 12 u xxz + ^“xyz^ + 7 ^22 u xxxz + 2I 13 u xxyz + ^“xyyz^ 


u(R) _ t u 

V y J {I 02 u v 


V (R) -£ fi u 
v z " ^ (I 02 u zz 


( r ^)[Ii2 u X yz + ^“yyz^ + 7 ^22 u xxyz 4 2I 13 u xyyz + ^Vyyz^ 

(E.9) 

2 

( 7~ )[I 12 u xzz + ^“yzz^ 4 7 ^22 u xxzz 4 2I T3 u xyzz 4 *04 u yyzz^ 


( r Q ) [I 30 u xxz 4 ^iScyz 1 4 2 ( r Q ) [I 40 u xxxz 4 2I 31 u xxyz 4 ^“xyyz^ 


v ( p ) . t 4 n 
v y ' ~J {I 20 U ) 


V z P) = \ {I 20 u zz 


" ( F^ [I 30 u xyz + ! 21 u yyz^ + 7 ( H [I 40 u xxyz + ZI 31 u xyyz 4 J 22 u yyyz^ 
° ° (E.10) 

t It 2 

if-) t^o^zz 4 ^“yzz^ 4 7 ^40 u xxzz 4 2I 3T u xyzz 4 I 22 u yyzz^ 


“7 {I 10 u x 
r o 


t 1 t ^ 

)[I 20 u xx 4 ^^xy^ 4 7 ^30 u xxx 4 2I 21 u xxy 4 I 12 u xyy^ 
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v^ } - ^ {I 10 u y - (^~) E 1 20 u xy + ! n u yy ] + 7 [I 30 u xxy + 2I 21 u xyy + l ^m ]) 
r ° ° (E.ll) 

v‘ ,x) - 4 dlo u z - ‘7- )[I 20 u xz * 'llV 1 * i [I 3° U X» 2 * n 21 a xyz 4 
^ r o u 

r o 

v‘ ,y) ■ 4 <‘<n“x - (f-'C'n-xx 4 >02“xyl * i t ,2[I 2’ u xxx + 2Il2 “xxy 4 'oAw 11 


»ny) . 4 <i „ - 4 I 0 2 u yy ] 4 * '^'“xxy * 2I l 2 u xyy 4 •Wniy 1 

r o ° (E.12) 

2 

,(>yl . 4 (i 01 » z - Ir>"n“K! 4 I 02 V 1 + 2 ( 4 ’ tl2, “xx* 4 2I i 2 “«yz 4 ! 03 u yyz : 
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APPENDIX F 

FAR-FIELD SOURCE FORMULAS 


As usual ii denotes the unit normal vector to a projected flat panel and T E 
and je are, respectively, unit vectors along the x- and y-axis of the panel 
coordinate system (Appendix B). Define r Q as the vector from the origin of 
panel coordinates to the point where velocity is being evaluated and r Q as 
its magnitude. Certain auxiliary vectors are needed: 





The far-field expressions for the source velocities are 


tfO) 

f(P) 

*<R) 


t^I r 
l *00 r o 

2 r 
r o 


t 4 

r o 


= -"3 l U D 
r o 


t 4 * 

“J *02* 
r o 


(F.l) 


(F.2) 


(F.3) 
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(F.4) 


fl(lx) 


t r n 

T — 

10 r o 
r o 0 


% [ Vi 


* >nty 


jdy) = 


, ?. t « 


[ ! 1 + *02 V 


The fact that three auxiliary vectors (F.l) are required means that in effect 
a transformation into panel coordinates has been performed. It appears to be 
somewhat faster computationally to use the vector far-field forms above rather 
than simply truncate the intermediate field formulas in panel coordinates. 
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APPENDIX G 

SOME SPECIAL NEAR-FIELD FORMULAS 

The near-field formulas of Appendices D and I can have numerical difficulties 
under certain circumstances. Some of these are inherent in the formulas, while 
others are due to extremes in panel dimensions - particularly very long thin 
panels. Special formulas have been developed to deal with these situations. 

"Small Logs" (Long Thin Panels) 


First consider side 32. If 


r 3 * ' fe < c, ■ ID' 3 

r 3 + r 2 * d 32 ^ 


(G.l) 


Then define 


.2 _ "32 A 

E '^T 


d 3 = 


— (x - e 3 ) + - 1 _ (y - V 


^ + 4< 


J\ + 


d 2 = d 32 ' d 3 


In the argument of the logarithm L^2) set 


1.1 1.2 1 . 1 __ 1 » 4 

r 3 * r 2 ' d 32 * 2 % * ' 8 ( „3 * d | )e 


(G.2) 


(G.3) 


For Side 41 

A similar procedure is used for r^ + r^ - d^. Define 

h 2 

2 _ h 41 + z 2 

E T + z 


1 + m. 


41 


(G.4) 
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d 4 “ 


m 


41 


A77 


(x - ? 4 ) + 


'41 




(y - 1 ) 3 ) 


(G.4) 


Then if 


d l = d 41 " d 4 


r 4 * r l ~ d 41 , £ 
r 4 * r l + c vT 3 


(6.5) 


set 

r 4 * r, - d 41 ^ - J (Jj 

The formulas are slightly different for the other two sides. If 

r l+ r 2- d 12 <£ 
r ! + r 2 + (* 12 3 


(G.6) 


(G.7) 


define 


e 1 2 = (y - n ^ 2 + z 2 


d l = x " h 


(G.8) 


d 2 ' d 12 " d l 


and set 


1 .1 1.2 1 .1 1_. 4 

= 9 (— + — ) £ " 8 .3 e 

1 G 2 


Vz- d i? ‘ 2 T d 2 o d . „ 


(G.9) 


Finally, if 


r 3 * r 4 - d 34 . . 
r 3 * r 4 + d 34 3 


(G.10) 


define 


e 2 = (y - n ,) 2 + z 2 


d 4 = (x - r 4 ) 


(G.ll) 
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d 3 = *34 " d 4 


and set 


+ r, - d 


34 


1 

1 


,1 A 1 % 2 
( ^ + ^ )e ' 


/l . 1 % 4 
( 3 + 3 )e 
d 3 d 4 


(G.12) 


"T-Derivatives" 

The individual T-derivatives in Appendix D become indeterminate if the point 
(x,y,z) is on the extension of a side of the panel. The "fix" of Ref. 5 is 
designed to remedy this, but it is inadequate for two reasons if the panels 
are long enough. One is that if the point is near the side itself, the fix is 
inappropriate. The other is that the criterion is too stringent if the point 
is indeed near the side extension. The following appears to be a reasonable 
-fix of the fix." 


1. Slanted sides 


If (y - hjMy - n 3 ) < o, skip this part. Otherwise calculate 


2 h 2 . ,2 

e 32 = h 32 z 


'41 * h 41 + z2 


(G.13) 


If e| 2 /y 2 < °* 00 °1, use 


3T< 32) + 3T< 32) 


3x 


3T* 32) 3T 3 

5z + Tz 


dx 

(32) 


3T< 32) 3T< 32) 

2 + — rrt — = 0 


3y 

"*32 


3y 

i 


A 


+ m: 


I y - n 3 1 " ly - n 1 1 


32 


If e^/y 2 < 0.0001, use 


(G.14) 


3T 


(41) 


3T 


(41) 


3T 


( 41 ) 


3T 


( 41 ) 

4 


ST 


ST 


ST 


+ z_ * 0 


(G.15) 
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9T ; 41) 9 tJ 41) 

— 3 T“ + — XT' 


m 41 i 1 _ 1 

l — r ly ‘ V ly ' "3 r 

/' + "41 


Parallel sides 

If h 32 h 4 1 < °» skip this part * 0therwise calculate 


? 2 2 
e ]2 = (h - r^) + z 


e 34 - <y - ^ )2 + z2 


(G.16) 


If 


2 

e 12 


{x - + y/2) 


- 2 < 0 . 0001 , 


(G. 17) 


use 


3T< 32 > al< 4, > 

ax““ + ax 

ai< 32 > + 3 T< 4T) 


(32) aT (41) 


air " 7 aT 


1 


ay ay 


= 0 


"41 


w 32 


az az 


l x - YfT ' lx - s 2 l 


(G. 18) 


If 


2 

e 34 


{x - % + V/2} 


,< 0.0001, 


(G.19) 


use 


3T< 32 > + aT< 4 " 


3x 


8 x 


aT* 32) aT^ 41 ) 


9 T< 32 ' ar< 41) 


az 


1 

~5F 


ay ay 

•"32 m 41 

-nr^-qT‘Tx^T 


{G .20) 
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Edge Vortex Formulas Near Extended Line 


These are needed for small values of 


q 2 = (y - iij ) 2 + z 2 


(G.21 ) 


If 


J 7- < 2.5 • 10' 2 


' ^1 or 2-* 


(G.22 ) 


use the following formula for Jg n (F) 

sgn(x - £,) , , 

•V F > * ri-T — «• 7 TPT * ~ 71 PT 1 


{^2 ~ x ) i 


Ui * x) 


n(n - 1 ) „2 i 1 

" 2 ( n + IT q L 


1 


( C 2 - x ) 5TT ' u , - x)" +T Jl 


(G.23 ) 


"Simpson’s Rule" 

Define R as the distance of (x,y,z) from the closest point of the line vortex 
and d as the length of the line vortex. Then if 

R/d > 10 (G.24) 

do not use the recursion formulas of Appendix I, but calculate the J mn by 
the three-point Simpson's rule 

5 , ^ - h h * l 2 

J mn = J 2 FU)d£ = 2 - 6 -~ 1 LFUj) + 4F ( — 2 -~) + F(^)J (G.25) 

^1 
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APPENDIX H 

CALCULATION OF VORTICITY INDUCED VELOCITIES IN TERMS OF 
SOURCE INDUCED VELOCITIES 

The calculation of the vort icily influences can be made much more efficient by 
expressing them in terms of the corresponding source influence, which of 
course, must be calculated In any event. The use of this procedure was put 
forward in Ref. 21. The portion of the theory that is needed for the present 

purpose is quite easy to state. 


Suppose there is a variable source density o on a portion of a plane or 
curved surface S. The velocity due to this at a point (x,y,z) is 

^ (source) = H cdS (H,1) 

where ? and r have their usual meanings. If there is a vorticity distribution 
on S of strength 

5- a* 

The Biot-Savart law gives the resulting induced velocity as 

? (vorticity) = JJ ^ (H,3) 

S r 

Then if t is a constant vector and if w has the same spatial variation as 
o, the velocity due to the vorticity distribution may be expressed in terms 
of the velocity due to the source distribution as 


1 (vorticity) = t x (source) (H.4) 

since £ can be resolved into components, each of which has a constant 
direction, the restriction to a constant £ is not serious. Although the above 
results apply to a curved surface S, it is far simpler to apply to a flat 
surface. In the present context the above is applied to the flat projected 

panel. 
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Figure 8 illustrates the projection of a curved panel S on the surface of a 
flat panel A in the tangent plane. In particular. Figure 8 illustrates r from 
a point of S and the vector r f from a point of the projected flat panel to 
the point (x,y,z). Evidently 

r f = (x - E)t e + (y - n)j e + z£ e (H. 5 ) 


where I . j , lc are unit vectors along the axes of the panel coordinate 
e 7 *e * g 

system. As in Appendix A, the vertical distance s between the curved panel 
and its projection is approximated by its leading term which represents a 
surface of second degree 

C 2 = PE 2 + 2QEn + Rn 2 (H.6) 

The aim is to obtain a consistent two-term expansion of Eq. (H.3) and express 
the results in terms of source effects. From Eq. (2.6.7) it is seen that a 
two- term expansion of the vector vorticity distribution is 


U> = U) + 


“I 


where 


is zero order and 


U) 


= y 1 

y e 


Ve 


(H.7) 


(H .8) 


. 9 ) 

is first order. The constants P x »u xx » etc. are the derivatives of the equiv- 
alent dipole distribution as given by Eq. (2.6.4). From Fig. 8 it can be seen that 

r = r r - cJc (H.10) 

i c “ 

Thus a two-term expansion ot the vel' Hy at (x.y.z) due to the vorticity on 
the panel is 


U)i 


- 2 <V + 


v yy n)i e 


- 2(Vx X E + vixynjjg + 2[-(QE + Rn)^ + (PE 


Qn)vy]k 

(H 




I 1 


u) x r 


dS = 


i 1 


w X r- 
r O f 


(1 + 3z — w ) + 


-m- -r t jr 

w. X r, w x k 

- L X 1 -^- 9 TT 1 ^ dA 


(H.n ) 
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- , ( r f f n + 3 Z i> - sAw ! /- 

1 7' s 


X r. 


dA 


. , *. ho x(0) an( i terms of the source 

By taking the graven of the ♦ x n can be seCT that the 

expansion in Eqs ^ u just the sun, of superscript 0 and c 

ZZ\ e Z for unit source density. Specific,,!, this cognation is the 

velocity 


f. = flO) a [P» (P) t 2Qf <Q) ♦ R» lR, l 


(H.12) 


rn« (241) To analyze the last term of 
This same combination appears in Eqs. (Z.4.U. 

Eq. (H.ll), collect terms in Eq. (H.9) to obtain 


w 1 * ZUi x + ^ y ) 


(H. 13) 


where the vectors 


\ = ^xy’e ' v xx J e + 


(Pv - QV*e 


Oy • Ve - Ve * «S • '“A 


(H.14) 


- —• » - rr: “ 

. «— >•-. -* •- 

;ux) and ;ny) Of Eq. (2.4.D- 

nos the velocity due to vorticity on the pane, «, be expressed in terns of 
source velocities as follows 


j =;„<«■ 

U) O 


♦ 2[q 


x V (U) ♦ q y X 


^ (ly) ] 


(H.15) 


. * that Eo (H 15' can be evaluated directly in refer- 

It is interesting t. «W ttat E, < , hm calC u,.ted and 

TT hit sysur ith r^d to tne velocities due to the vorticity. this 
put into this system, *n y reference coordinates 

r. - 

neve r arise. If the source velocities have been computed by far 
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formulas, they simply are used in Eq. (H.15), so that in effect the vorticity 
calculation uses the source far-field procedure. The present code takes 
advantage of the second of these facts, the use of far-field source formulas, 
but performs the calculation in panel coordinates. 

The implementation in panel coordinates proceeds as follows. 


Define 


“1 = “ vT 


“2 = vT 


hr- 

*1 " vT * c{ri l + V 
h s 

e 2 = - vr + c(nj + *i 3 ) 


then 


B F (i e * J e t-v-e,] ► k e [#,»; * «,*;]) 

+ B s {t e [-»V 2 ] * 3 e t-V-B 2 ] ♦ t e [6 2 v- + a 2 v;]) 


Define 


6 1 = 2w 


6 2 = " 2w 


e l = c 


ei =-c 


then 


^x x 


(H.16) 


(H.17) 


(H . 1 8 ) 


* (lx) = Bp {T e [Vy 1x) (PB, - Qo-j ) + J e [V i 1X)(P6 l * ^ " V z' X ' 6 l ] 


i(lx) , 


+ Ic e [V^ lx) 6 1 ]} 


y 

• e [v y 

"S ' 'e [V y 

lx), 

e 1 - y C 2 


+ Br ft [vJ lx) (PP 2 - Qa 2 ) + J e [V^ lx) (PB 2 - 0«2 > - 


r(lx). 


(H - 19 ) 


+ IT rvi lx) 6,]} 
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q y x tf (1y) * B p {T e [-V^ ly, 6 1 - V y ly) (QB-j - Ra-j) + J e [vJ ly) (QB 1 - Ra-j ) - V* 1 ^] 

+ k> (iy) £i + v (iy) 6i]} 

(H.20) 

+ B S {T e t- v z ly)s 2 ‘ v l 1y)(QB 2 ' Ra 2 } + V V x ly)(QB 2 " Ra 2 ) " * 2 ] 


♦ ^ V i ly)£ 2 + V i ly)6 2^ 
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APPENDIX I 

FORMULAS FOR THE EFFECT OF NEAR-FIELD LINE VORTEX 
ALONG A STREAMWISE EDGE OF A PANEL 


Derivation of the Influence of an Edge Vortex 

The equation of the curved panel Is Eq. <H.6>. For definiteness consider the 
case when the edge in question lies in the plane n * n,, the f1rst 
M-line (Fig. 3). The modifications for the case of the second N-line are 
obvious. Thus the curve c along which the vortex lies is 

C = c(U • P? 2 + 2<Kn, + *n 2 lI,1) 

The unit vector along this curve is 

[1 ♦ Tl e l »- 2 > 

/ 1 + T 2 


where 


T = 2(P£ + Qni) 


(1.3) 


The velocity due to the vortex is 

? = J £ *il ydS u* 4 > 

r c r 3 

where y is the edge value of the equivalent dipole strength. Arc length 
along the curve is related to distance in the tangent plane by 


ds = / 1 + T 2 dS (I * 5) 

Thus with r expressed in panel coordinates (Fig. 8) 

(t x r )ds * (C 0 -Tty-t^) ]* e 

[ z + T(x - E) + c]J e (I * 6 ' 

[(y - tvj) + 0 
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where the terns In the first coluw of Eq. (1.6) are first order, and those In 
the second column are second order. This expression is exact except for the 

approximation x, = c 2 . 


3 , 


As shown in Appendix A, a three- term expansion of 1/r is 


[1 + 3c 1 + 3(c* + c 2 )] 


(1.7) 


where r f is distance from a point of the flat tangent panel and where 

c 

c i ' zr 


3 2 1 *2 

c 2 = i c r?i 

r f 

Along the N-line the equivalent dipole strength varies linearly 

y * B(h + K) 


( 1 . 8 ) 


(1.9) 


where h is the total arc length along the N-line up to the n-axis of panel 
coordinates (see Section 9.2 of Ref. 2), and B is the unknown value of vortic- 
ity that is determined from the Kutta condition. The fundamental flow is 
obtained by setting B equal to unity in Eq. (1.9). Multiplying the above 
expansions gives the components of the vortex velocity as follows 

V s JL {o - [T(y - iVj)h] - [T(y - n-|)(3c^h + £)]}d£ 
r x 5, r 3 


y = J ^ {-zh + [-z(3c.jh + E) + h(T(x-j— S)i+ 

r y h r f 


(I. 10) 


+ [-z{3(c2 + c 2 )h + 3c,*} ♦ C3c,h + E)(T(x - l) + 

V = (y - h.) r 2 \ & + t 3c i h + ^ + [3{c l + c ? )h + 3c l*’ ]}dr ’ 
r z ^ r f 

The integrals in Eq. (I. 10) have the form 
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(1.11) 


mn 


i ' 2 i « 

i r f 


Once the J n „ and J ln have been calculated the others are calculated from 


'On a,,u ''In 
the recursion formulas 


J mn = J (m-2)(n-2) + 2xJ 'm-l)n ' p J 


(m-2)n 


( 1 . 12 ) 


where 


p 2 = x Z + (y - tv,) 2 + z 2 


(1.13) 


The required Jg n and J-j n are 


. . r l +r 2 td 12_ ,(12) 

J 01 * 109 r, 4 r z - d" ' L 

J„ • r 2 - r, ♦ xJ 0 , 


1 r 2 


J 03 * ^ r, 

q L 

J 13 = rj ' r“ + xJ 03 


U - x - x 


(1.14) 


5, - x C, - x 
Jnc = J ^ . 

3q Z r 2 


1 r ” 2 ~ '1— +2J 03 ] 


T 

r 1 


J 15 * " 1 7 ] + XJ 05 


r 2 r l 


,07 = i? 


1 r s 2 




+ W 05^ 


1 ,1 1 


J 17 - * 5 <3 ' 3 ’ + ^7 


r 2 r l 


where 


q 2 = (y - Oj) 2 + z 2 


(I.T5) 
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and where r } and r £ are, respectively, distances of the point (x,y,z) from 
the ends of the interval, i.e., 

r 2 - (x - 5 k > 2 * (y - V 2 * z 2 »•«> 

In terns of certain auxiliary functions, F n , the velocity components of Eq. 

( 1 . 10 ) are 

V rx = " (y " n l )[hF l + F 5^ 

V ry * -zhJ 0 3 ♦ [-ZF 2 * hFj] - z[3hF 4 + F 6 ] + F 7 U.W) 

V rx = + * y “ n l^ hJ 03 + F 2 + 3hF 4 + F 6^ 

The auxiliary functions for on-body panels are: 

F i * ZPJ'jj + 2 Qn-|JQ 3 

F 2 = 3zh[PJ 25 + 2Qy 15 + R^J q5 ] + (J 13 > 

_ 2 , 

F 3 = 2PxJ 13 + 2Qn 1 xJ 03 - P0 Z 3 + Rn l°03 
F 4 = I z2q 7 " \ Q 5 

Qj = [P 2 J 4j + 4PQn 1 J 3j + (2PR + 4Q 2 )ni0 2 j + ^ RT 1 l J lj + (1.18) 

F 5 = 6 zh[P 2 J 35 + 3PQn-j J 25 + PRn 2 J 15 + 2Q 2 t^J 15 + Q R 4W + {2PJ 23 + 2Qtv I J 13 > 

2 

F g = (3z[PJ 35 + 2Qn-,J 25 + Rn 1 J 15 ]> 

F ? * 3 zh[-P 2 J 45 + (2P 2 x - 2PQn 1 )0 35 + 6 PQxn 1 J 25 + (4Q 2 xn 2 + 2 PRx 2 + 2QRn 1 )J 15 

+ (2QRxn 2 + R 2 n^)J 05 l + {“F'ljj + ^Px^S + (2Qxn 7 + Rn^ )■! 13 ) 

The formulas for wake panels are obtained from Eq. (1.18) be deleting all terns 
in {} and replacing h by L (total). 
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d 


For the semi- Infinite last wake additional changes are made to the formulas 
(1.14) for the J mn corresponding to 

Z 2 - • r 2 -*• ® S 2 /r 2 - 1 (1.19) 

Furthermore, P and Q are set equal to zero. 
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APPENDIX J 

FAR -FIELD EDGE VORTEX FORMULAS 


Compute 


rp • [x 


C-. + C- ^ ? 2 

( 1 r — )] + [y - + z 


(J.l) 


If 


use 


where 


u 2 - t,) 2 /r* < 0.001, 

V rx = -<y - n,)T 0 I 

5, * E, 

v ry = t_z * T o (x ' ' 2 * 5 o ]I 

v rz = {y ‘ V 1 




T„ * 2[P 

Cl + c, 2 


Ci + ^2 
1 g - + Qn-j 3 


c 1 + c 2 


C 0 = (P — 5 — ) + 2( *V 7 


) + Rni 


( J.3) 


For the second N-line the obvious quantities are replaced by the corresponding 
ones. The above equations replace the elaborate formulas of Appendix I. 
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APPENDIX K 

PARABOLIC CHORDWISE VORTICITY 


The assumption of constant .orticity around a wing section (linear .ariation 
of the underlying dipole strength) can lead to numerical difficulties i on cer- 
tain wings with very thin trailing edges. Accordingly, a second hordv se 
variation option is required. The important consideration is to have 
•bound- vorticity strength approach zero at the trailing edge on both upper and 
lower wing surfaces. This is accomplished by a quadratic global variation of 
vorticity 5 as a function of arc length along an N-line. While only two , local 
chordwise variations have been incorporated into the present method many such 
variations are possible. As will be seen below, the required modifications 
the program are quite minor. This flexibility is due to the use of vorticity 
as an auxiliary singularity. 

To implement the parabolic vorticity option, the linear variation of the under- 
lying dipole strength along an N-line is replaced by a cubic variation having 
zero derivative at the upper and lower trailing edge. Specifically, 


L . 

v • Bs < 3 l (total r ' 2 '■L (total ) 1 1 


(K.l) 


where s is arc length along the N-line. The above is a global variation. The 
variation over an individual panel can be no higher a degree than quadratic, 
and in the present method has been taken as linear. It is assumed that the 
underlying dipole distribution on a panel agrees with Eq. (K.l at the corners 
of the panel and varies linearly in between. Thus, the overall behavior 
that of an inscribed-polygon approximation to Eq. (K.l). 

For an individual panel the arc length measured along an N-line from the 
trailing edge to the lower corner of the panel is ly + E, or h $ + E, while 
arc length associated with the upper corner is h F f Eg or h $ + Eg. The linear 
function that agrees with Eq. (K.l) at these two values of arc length is 


# \ 


n / ii 
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where on the two N-lines the constants H and I have the values 


H 


3 (hr ' 5,5,) - r .. g . '4 ' ^l^Z (3h F + 5 1 + 5 2 )] 


F - L f (totaTT '"F ' M'2- ^ , tota „] 


(Zh r - ?,+ l 0 ) - 


2 [3h p + 3h p (^ + K 2 ) 


h = Lp (totiTT '‘"F ■ T '2' [Lf (total )] 

+ A*4* 5,5 2 J 

H S - rp^U "4 - - 1 , t ou T? [ 4 - ( * t3h s * *» * 


(K.3) 


Ic = 


(2h c - Z.) - 


S = L s (totaTT l4,, S - ^3 -4' [L {total)] 


T [3h| + 3h $ (C 3 + C 4 ) 


♦ 4 * 4 ♦ V4 ] 

where all symbols have the sane Meaning as in Section 2.6.1. Thus the 
quadratic form (2.6.9) for the variation of dipole strength over a panel is 

replaced by 

Bplp - B S I S 8 F Hp - B S H S Bslgn, - Bplp^ BsVl - BpH F n 3 
u = Cn + # n + w ^ w 

+ c(B p - B s )(n - n 3 )(n - Oj) (K * 4 

The dipole derivative formulas of Section 2.6.1 are modified in an obvious 
way, specifically 

a. Terms containing c are not changed. 

b. In terms containing h p or h $ these quantities are replaced by 
H p or H s , respectively. 

c. All other terms are multiplied by I p or I s as appropriate. 
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The wake formulas are unchanged except for c. In the constant chordwise 
vorticity option, the parameter c is nonzero on wake panels if the "piecewise 
linear" spanwise variation of vorticity is used. However, if the parabolic 
chordwise vorticity option is used, c is taken as zero on all panels. 

The near-field edge-vortex formulas (Appendix I) are modified as follows: 

= unchanged 

F 2 = replace h p by H p 
multiply J 13 by I p 

F 3 = unchanged 

F 4 = unchanged (k.5) 

Fg = Replace h p by H p 

multiply [2PJ 23 + 2Qn T J 13 ] by I p 

Fg = mult ply entire term by I p 

F 7 = replace h p by H p 

multiply [-PJ 33 + 2 P x J 23 + (2Q x n 7 + Rt^)J ] 3 ] by I p 

The wake formulas are unchanged. Replace H by L p (tot.) and set I = 0. 

Note that the terms multiplied by I are exactly the terms neglected in the 
wake. 


In the far field (Appendix J), the only change is that I becomes 

dp S] + S n 

1 = “3 (H F + 2 V 

r F 
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APPENDIX L 

B DERIVATIVES AT SECTION EDGES 


If the k-th strip is last in a section, the square bracket in Eq. (2.6.30) of 
Section 2.6.5 is replaced by 


[DB k _ 2 ♦ EB^ ♦ FB k ] 


(L.l) 


w,. 


D = 


L k “ 1 ! W . k -) 


Vl + l/zl V2 + \ ) \-2 +W k-l 


(L.2) 


r . 4w \-i* i/2 ;v2*v 
k (w k . 2 + Vl M Vl + V 

V 2 + 3 V1 + 2w k 
F = w k ( Vl 7 v L Vi + 

If the k-th strip is first in a section, the square bracket is replaced by 


[DB k + EB k+1 + FB k+2 ] 


(L.3) 


D = -w 


2>i k + H+l + "k+2 


k l« k + Vl’ l Vl ♦ '^'“k + V2 I] 

V, *m(\* y 2 > 

' k >“k * Vi"\.i + V2’ 


(L.4) 


w k ”k * Vl 

r " * ' Vl + 1/2 K * V2 1 ‘ Vl * V2 


If a sec Jon has only one strip, eliminate the square bracket, i.e. 


D = E = F = 0 


(L.5) 


If the section has two strios, use 


D = 0 


L-l 
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E = - 


(L.6) 


2w k 

w k+l + w k 


F = -E 


for the first strip, and 




+ Vi 


for the second strip. 


(L.7) 
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APPENDIX M 

CONVERGENCE ACCELERATION SCHEME 


After each iteration, a convergence acceleration procedure is invoked in which 
a new solution is defined in terms of a linear combination of the previous 
solutions. This appendix will give the details of the calculation of the 
required linear combination. 

Let X°, X 1 , ... X k and RES 0 , RES 1 , ... RES k be k + 1 successive approximations, 
and their corresponding residual vectors. Since there are N+L unknowns to be 
solved for, we define the (N+L) x (k+1) matrix whose columns are the solution 
vectors 

[X] = [X°, X 1 , ... X k ] (M.D 

Similarly, define the residual matrix 

[RES] = [RES 0 , RES 1 , RES 2 , ... RES k ] (M-2) 


Define the "linear combination vector" 

f 


such that a new approximation X' is given by the matrix product 

f ol 

! 

X' = [X] 


and the corresponding residual is 


f o 


RES' = [RES] . 

if. 


(M.3) 


(M.4) 


ri-i 
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Note that we can write 
RES’ = R - AX' 


= R - A[X] 


= R(1 - I fj) + t(R “ AX°),(R -AX ), ... (R-AX )] 
i-0 1 


k 

i_ j 


= R(1 - I f,) + [RES 0 , RES 1 , ... RES k ] 
i=0 


Lv 


Therefore, by choosing 


k-1 




i=0 


the first term will disappear, and a new residual is given by 

f 


RES' = [RES] 


k-1 

1 " f o ‘ *** ” f k-l 


Define (k+1) x (k+1) "modified- unit matrix I] by 


(H.5) 


(M.6) 


(M.7) 
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I 

1 


•1 -1 -1 .. 


-1 1 


from which the linear combination vector can be written 


f- -1 


*0 



f »1 


• 


F 

• 

■ >1 

• 

• I 

* X 1 

1 

• 



_ f k_ 


f k-l 


- -* 



_1 




where F is the (k x 1) vector 

“ f o ' 
F = 


The new residual vector RES' can now be written 
the old residual matrix and the unknown vector F 


as a matrix product 
in the form (using 


RES' = [RESHl-j] 



Define the norm of this vector, I IRES || by 


I IRES' 1 1 2 


T 

RES' • RES' 


= [F T 1] [I,1 T tRE$f [RES] 


Nil 
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(M.8) 


(M.9) 


involving 
Eq. (M.4) ) 

(H. 10) 

(M.ll) 




The right-hand side of this equation is a quadratic non-negative (scalar) 
function of the k unknowns f Q , f^ ... f^.-j. The minimum value must therefore 
occur at the point at which 

•af: I IRES' 1 1 2 ■ 0 for i = 0, ... k-1. 

i 

This will therefore provide k linear equations which can be solved to minimize 
Eq. (H. 11) . 


Calculation of Partial Derivatives of 1 1 RES * 1 1 2 


First define the sywnetric matrix P by 

P = [RES] 7 [RES] 

i.e. P = [P.j], where P.. = RES 1 • RES J (M.12) 


(scalar product of 1*^ and residual vectors). Now partition the matrix P 
in the following manner: 


P 11 

P 12 

(kxk) 

(kxl) 

P T 

i2 

P 22 

(lxk) 

(1x1) 


(M. 13) 


so that P^ = scalar products between all of the residuals RES 0 ... RES k_1 , 
while P 12 consists of scalar products between latest residual RES k , and aJ2 of 
earlier residuals, RES 0 ... RES k '\ and ? 22 is scalar ||RES k || 2 . 


Next define the symmetric matrix Q by 


Q - [I^ 1 [P] [I'], (H. 14) 

and again partition the matrix Q to separate the last row and column: 
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I 

Q 11 

— 1 

*12 



-1 


1 

(kxk) 

(kxl ) 


1 I 

• 

• 

• 


Q = 1 

<>; 2 

Q 22 

_ 


-1 



(Ixk) 

(lxl) 


[0...0] 

1 


Straightforward matrix multiplication c 


n 


12 


12 


22 


(M.15) 


[- 1 ...- 1 ] 1 


or 


and 


or 


and 


Q n = ? u + P 12 H1 T + [ ' 1]P 12 + p 22 [ ' 1][ " 1]T 

q ’V Pi v Pl2 i' p V P22 


(M. 16) 


Q 12 = P 12 + P 22 [ " i:i 

q 12 1 = " P 22 

^22 = P 22 


(M.17) 


Eq. (H- 11) nw can be written 

| |RES' 1 1 2 - [F T I) 


<11 




*12 


■‘22 


L 


(M. 18) 


,„d partial differentiation of this expression with respect to f„. V" 
f provides a set of k linear equations: 

Qn F + Q 12 = 0 

„ - (Me 19) 

or Q 11 F * ' Q 1Z 

the solution of which provides the unhook k acceleration coefficients. 

Given this solution, we can define the acceleration vector 
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from which the new solution X' is given by 

X' = [ X ]F * 
and 

RES' = [RES]F‘ 


(M.ZO) 


(M.21) 
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appendix n 

the compressibility correction 


The user inputs the average ^compressible velocity V„ -H* 
pond to physical conditions in the region ’ 
equal to ai and the average density ratio e ,s 


e = P_ = 0.6339 

Pt 


(N.l ) 


„ this occurs the point is labeled -choked- on the output. 

„ Vl < ,i. it is used as it stands to co^ute (o/P t ) » « 
procedure. The iterative equation is 

v 2 2.5 

1 /i , I -i (N.2) 

•-i W V 7 1 

with initial value e * 1- 

Finally, the co^ressible velocity magnitude V is calculated free, 

tf ■ 


1 •" 

V - Vi (i> . 


v i 

“X 


(N.3) 


V 1 is the magnitude o. the local equivalent incompressible velocity 
where V- , woinritv is not changed. 

(Section 3.3). The direction of local velocity is n 

. n c , n that an equivalent incompressible 

TveX ZZ't'ZZS* - - control station, for con^atibility at 
the control station, the input should insure that 

(N.4) 


V = V' 
v i c 


or the computed surface ^ ”,^1, contradictory 

— *• * the re9lon of interest ,s 

from the control station.. 
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APPENDIX 0 

OPTIONS OF THE COMBINATION PROGRAM 


0.1 Incompressible Option 

If the flow is incompressible, this option is selected and only the following 
quantities are input: 

V freestream velocity 

co 

V average velocity at the control station 

c 

V , reference velocity used in computing pressure coefficient 
ref 

a angle of attack 

B angle of yaw 

0.2 Freestream Conditions 

For compressible flow the freestream conditions are defined by inputting angle 
of attack a, angle of yaw B, and three additional quantities: 

either velocity or Mach number 

either total pressure P t or static pressure P $ 

either total temperature T t or static temperature T $ 

Then the preliminary calculations are as follows: 

0.2.1 M ot Input 

(a) If P t is input, P $ is given by 

i ? -3/5 

p s sP t ( 1 + l , £ ) (0J) 
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( 0 . 2 ) 


If P s is input, is given by 

p t * P s" + 5 H ~> 

If neither is input, the default is 


3/5 


and P s is as above. 


P t = 2116.23 


(0.3) 


(b) If T t is input, T s is given by 

-1 

T S * T t< ,+ 5>£> 


(0.4) 

If T,. is input, T t is given by 

T t = T s ( l + ^ M i) 

(0.5) 

If neither is input, the default is 


T t = 518.67 

(0.6) 


and T s is as above. 

In either case stagnation and freestream sound speeds a^ and a s are calculated 
from 


a t * 49 A t , 


*$ * 49 *1 


and V ro from 


V = a.M (1 + i 

co t ® 0 00 


- 1/2 


0.2.2 Input 

(a) If T t is input, a t is given by 


H 5 49 A t 


Ho, is then calculated from 

V 


l v m 2 


- 1/2 


(0.7) 


( 0 . 8 ) 


(0.9) 


( 0 . 10 ) 
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The remainder of the calculation proceeds as in 0.2.1 above, 
(b) If T s is input, a s is given by 

a s - 49 /T s 

and Mod by 

M = Y/a, 

co oo 5 

The remainder of the calculation proceeds as in 0.2.1 above 
0.2.3 Additional Freestream Quantities 
Built into the program are constants 

g = 32.174 

R = 1715.63 

The following quantities are calculated: 


Total density: 


Static density: 

P s 

p s ' RT S 

Dynamic pressure: 

q. - °* 7p t { 

Pressure ratio: 

Ps/Pt 

Density ratio: 

Ps/Pt 

Temperature/sea-level ratio: 

. T t 

6 '5107 

Pressure/sea- level ratio: 

.. P t 
6 " 2116.21 


( 0 . 11 ) 

( 0 . 12 ) 


(0.13) 

(0.13) 
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Critical speed: 


a* ■ a t //T72 


Maximum velocity: ^max " ^ a t 

p s 

Equivalent incompressible freestream velocity: = V. (— ) 

Equivalent incompressible critical velocity: ai = 0.6339a*. 

0.2.4 Sunmary 

Three freestream conditons are input: or M^, P t or P s , T t or 

T s (or default values). Calculated and saved are 

V„. P t , P s , T t , T $ , a t , a s 
p t’ p s* q »* P s /P t» p s /p t’ 6 * 6 

I I 

a *» ^max* a * 

Nineteen quantities all together. 

0.3 Control Station Condition^ 

Input consists of one of the following three quantities: 

$ inlet mass flow rate 

V c average velocity 

M average Mach number 
c 

The remaining t» must be calculated plus some additional quantities. 
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0.3.1 V c Input 


p is given by 
c 


e c = p t 


1 V r 2 

P* [ 1 ' 3 ( "a^ ^ 


2/5 




by 


and M c by 


/y p OO 
« = gp c v c -ppr - 


V , v r 2 

u - !c [ 1 - i ( a 2 ) 3 

M c a t L » a t 


- 1/2 


0.3.2 M c input 
V c is given by 


» c - a t M c (1 + 5 "?> 

Then Pc and * are obtained as in 0.3.1 above. 


- 1/2 


0.3.3 * Input 


Here v c 


V r must be calculated 


iteratively by solving the equation 


Yi 75 

V c ” gUVpplkl/l^'Pt [1 ■ 1 /z(V c /a t ) 1 


starting with V c - 0. 

Once V c is known, M c and p c are 


obtained as in 0.3.1 above. 


(0.15) 


(0.16) 


(0.17) 


(0.18) 


(0.19) 
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0.3.4 Additional Control Station Quantities. 


These are calculated as follows: 


Pressure ratio: 
Density ratio: 

Dynamic pressure: 
Velocity ratio: 


V. 2 


3.5 


<W * [ 1 - F <a7> 1 


P c /p t 


"c ’ °- 7P t <fT> H c 


V«/V c 


Corrected mass flow: * cor = * -j- 

Equivalent incompressible average velocity: 


V c 


V (— 
c 'p t 


0.3.5 Sumnary 

One quantity, ft, V c , or M c is input. Quantities saved are 

V c* M c* p c* {P c /P t ] > {p c /p t ) * V ^ V a/ V c^’ "cor* V c 
a total of ten quantities. 


( 0 . 20 ) 


) 


( 0 . 21 ) 
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APPENDIX P 

ORGANIZATION OF THE INPUT POINT? 


The input to this program consists of the coordinates of a nunfcer of points. 
These points define the surface of the three-dimensional inlet around which the 
flow is to be computed. For the purpose of organizing these points for compu- 
tation, each point is assigned a pair of integers, m and n. These integers 
need not be input, but their use must be understood to insure the correctness 
of the input and to facilitate the interpretation of the output. 

For each point, n identifies the "column" of points to which it belongs, while 
m identifies Its position in the "column," i.e, the "row." The first point of 
a "column" always has m = 1. To insure that the program will compute outward 
norma! vectors, the following condition must be satisfied by the input points. 

If an observer is located in the flow and is oriented so that locally he sees 
points on the surface with m values increasing upward, he must also see n val- 
ues increasing toward the right. Examples of correct and incorrect input are 
shown in Fig. 34(a). In this figure the flowfield lies about the paper, while 
the interior of the body lies below the paper. Occasionally, it happens that 
despite all care a body is input incorrectly. If the entire body is input 
incorrectly - not some sections correctly and some incorrectly - the difficulty 
can be remedied by changing the sign of one coordinate of all the input points. 
This trick will give a correctly input body of the proper shape at perhaps a 
peculiar location. Otherwise, the input will have to be done over. If the 
inlet is input correctly (Step 2), but a cross-section (Step 4) is input so 
that its normal vector points upstream, the confcined flow will be correct, but 
the flux at the cross section will be negative. Clearly a control station with 
the wrong normal vector invalidates the calculation (Step 4). 

The body surface is imagined divided into sections, which may be actual physi- 
cal divisions or may be selected for convenience. A section is defined as 
consisting of a group of at least two n-lines. Within each section the n-lines 
are input in order to increasing n. On each n-line the points are input in 
order of increasing m. All n-lines in a section must have the same nunfcer of 
points, but this may vary from section to section. The first n-line of the 
first section is n = 1. From then on the n-lines may be thought of as numbered 



1 
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( ecutively through all sections, i.e., the numbering is r.oL begun over at 
the beginning of each section. Elements will be formed chat are associated 
with points on every n-line except those that are last in their respective 
sections. Points on these latter n-lines are used only to form elements 
associated with points on the next lowest n-lines. 

To illustrate this procedure, consider the plan view of a body shown in Fig. 
34(b). This body has been divided into four sections, as shown in the figure. 
The first section contains four n-lines, n = 1, 2, 3, 4; the second, five 
n-lines, n = 5, 6, 7, 8, 9; the third three n-lines, n = 10, 11, 12; and the 
fourth three n-lines; n = 13, 14, 15. The nuinber of points on each n-line are: 

Section =1234 
M = 4 7 4 2 

Notice that the line n =4 has only four points, the points m = 1, 2, 3, 4 and 
the m-grid of Section 1, which is listed in the figure along the n * 1 line. 

The lines n = 4 and n = 5 are physically identical. Some of the points on the 
two lines are physically identical but correspond to different values of m. 

This is of no consequence. In this scheme sections are completely Independent. 
No elements are computed corresponding to points on lines n = 4, 9, 12, 15. 

9 

There is no restriction that the m and n lines of different sections have to 
be roughly parallel. The arrangement shown in Fig. 34(c) is permissible. 
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